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Abstract 

Starting from an approximate microscopic model of a trapped Bose- 
condensed gas at finite temperatures, we derive an equation of motion for the 
condensate wavefunction and a quantum kinetic equation for the distribution 
function for the excited atoms. The kinetic equation is a generahzation of our 
earher work in that cohisions between the condensate and non-condensate 
(C12) are now included, in addition to collisions between the excited atoms as 
described by the Uehling-Uhlenbeck (C22) collision integral. The continuity 
equation for the local condensate density contains a source term which is 
related to the C12 collision term. If we assume that the C22 collision rate is 
sufficiently rapid to ensure that the non-condensate distribution function can 
be approximated by a local equilibrium Bose distribution, the kinetic equation 
can be used to derive hydrodynamic equations for the non-condensate. The 
ri2 source terms appearing in these equations play a key role in describing the 
equilibration of the local chemical potentials associated with the condensate 
and non-condensate components. We give a detailed study of these hydro- 
dynamic equations and show how the Landau two-fluid equations emerge in 
the frequency domain cjt^ <C 1, where is a characteristic relaxation time 
associated with C12 collisions. More generally, the lack of complete local equi- 
librium between the condensate and non-condensate is shown to give rise to 
a new relaxational mode which is associated with the exchange of atoms be- 
tween the two components. This new mode provides an additional source of 
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damping in the hydrodynamic regime. Our equations are consistent with the 
generaUzcd Kohn theorem for the eenter of mass motion of the trapped gas 
even in the presence of colhsions. Finally, we formulate a variational solution 
of the equations which provides a very convenient and physical way of esti- 
mating normal mode frequencies. In particular, we use relatively simple trial 
functions within this approach to work out some of the monopole, dipole and 
quadrupole oscillations for an isotropic trap. 



PACS numbers: 03.75.Fi, 05.30.Jp, 67.40.Db 
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I. INTRODUCTION 



At very low temperatures, the dynamics (collective oscillations) of a trapped Bose gas 
is described by the time- dependent Gross-Pitaevskii (GP) equation [|I| for the macroscopic 
wave function $(r,t) associated with the Bose condensate. At higher temperatures, when 
an appreciable fraction of atoms is excited out of the condensate, the dynamics of the 
trapped gas becomes much more complicated since it now involves the coupled motion of 
the condensate and the non-condensate degrees of freedom. In recent work ||^ the finite- 
temperature excitations of the condensate were calculated within the Hartree-Fock-Popov 
(HFP) approximation. This treatment of the excitations is equivalent to a description based 
on a generalized GP equation in which the condensate density ^^(r, t) is treated dynamically 
but the non-condensate density n{r, t) is treated as static. That is, one ignores the effect 
of non-condensate fluctuations on the dynamics of the condensate. Several authors |^,^ 
have recently extended this kind of HFP theory by treating the non-condensate fluctuations 
5n(r, t) using a time-dependent Hartree-Fock equation of motion (often referred to as the 
RPA). This describes the collective modes in the coUisionless mean- field regime. 

In contrast, in the low frequency collision-dominated regime, one expects a description in 
terms of hydrodynamic equations. In a recent paper , we discussed a method for obtain- 
ing two-fluid equations starting from a semi-classical kinetic equation for the distribution 
function /(p,r,t) of the excited atoms. The final result was a closed set of two-fluid hy- 
drodynamic equations in terms of the local densities and velocities of the condensate and 
non-condensate components. In this paper, we consider a generalized model and present 
a more complete discussion of the approximations which are involved. As we shall see, 
this new microscopic model, though approximate, captures all the essential physics of the 
collision-dominated hydrodynamic behavior involving the condensate and non-condensate 
components. 

Deriving two-fluid equations within the context of a specific microscopic model has several 
advantages: 

1. It allows one to see in a very clear fashion how the condensate and non-condensate 
degrees of freedom arise and how they are coupled. 

2. It gives a closed set of equations with all local thermodynamic properties and relaxation 
times determined in a self-consistent manner. It thus ensures that the static and time- 
dependent properties are treated in a unified manner. 

3. The equations we obtain can be used to derive the usual Landau equations (origi- 
nally developed to describe superfluid ^He 0,01) in a certain limit, and allow a clear 
identification of the Landau superfluid (normal fluid) density with the condensate 
(non-condensate) density. 

4. The description of the dynamics can be extended beyond the regime in which the 
Landau phenomenological equations are valid, i.e., when the condensate and non- 
condensate are not in local diffusive equilibrium. 
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5. An explicit microscopic formulation provides a starting point for improved descriptions, 
for example, the inclusion of the effect of the off-diagonal (or anomalous) self-energies 
which are important at very low temperatures 

In the ZGN paper f^, we only included collisions between the excited (non-condensate) 
atoms. Moreover, we only discussed the solutions of the kinetic equation for /(p, r, t) which 
correspond to dynamic local equilibrium within the non- condensate. This local equilibrium 
state is defined by the condition that the collision integral for the excited atoms (denoted by 
C22[/]) vanishes identically. As a result, no explicit dependence on the C22 collision integral 
appeared in the ZGN hydrodynamic equations. In the present paper, we make a major 
improvement by also considering the effect of collisions between atoms in the condensate 
with those in the non-condensate. This collision term (denoted by Culf]) allows us to 
discuss how the condensate is coupled to the non-condensate beyond mean-field effects, and 
how the condensate comes into dynamic local equilibrium with the non-condensate. In this 
connection, we should note that the different effects of the C12 and C22 collisions was first 
discussed in a qualitative way by Eckern |^ . The generalized set of hydrodynamic equations 
which arise when the effects of Ci2[/] are included will be referred to as the ZGN' equations 
to distinguish them from the set we derived previously . 

Even when collisions among non-condensate atoms are sufficiently rapid to enforce the 
distribution function / to take the form of a local equilibrium Bose distribution with chemical 
potential jl, it is possible that the non-condensate is not in complete local equilibrium with 
the condensate. In this situation, the chemical potential fl of the non-condensate can be 
different from that of the condensate, /ic. Thus the variable fidis = fj' — fi^c emerges as a 
key new local variable which is needed to represent the equilibration of the condensate and 
non-condensate components. In the original ZGN analysis, where the effects of C12 were 
neglected, the difference between /ic and /i was implicitly built into the ZGN hydrodynamic 
equations, but this fact was overlooked in our attempt to reduce these equations to the 
Landau two-fluid hydrodynamic equations. To arrive at this latter description, we assumed 
additional local equilibrium thermodynamic relations which implicitly forced fidis = (i.e., 
/i = /ic). Thus the ZGN equations are in fact distinct from the Landau two- fluid equations 
and, in general, will lead to different dynamical predictions. In the present paper, we show 
how the more general ZGN' equations reduce to either the ZGN or Landau descriptions, 
in the appropriate limits. This more complete description involves an additional equation 
of motion for fidis, which is coupled to the condensate and non-condensate velocity fields. 
The new relaxational mode which emerges involves the transfer of atoms between the two 
components, and allows one to understand the conditions under which the Landau two-fluid 
equations are vahd. 

The possibility of the condensate being out of equilibrium with the non-condensate was 
studied many years ago, mainly with regard to the behaviour of superfluid helium near the 
A point. In a well-known paper. Landau and Khalatnikov |p calculated the absorption of 
sound near a second-order phase transition due to the coupling to a relaxing order param- 
eter. Pitaevskii |T0|,[7| subsequently generalized the phenomenological two-fluid equations 



for superfluid helium to take into account the relaxation of the superfluid density near the 
A point. The same problem was looked at within the Mori formalism by Miyake and Ya- 
mada All of these approaches are phenomenological in nature since they introduce new 
dissipation parameters to describe the process of atomic exchange between the condensate 



4 



(superfluid) and non-condensate (normal fluid). In contrast, our theory is based on a well- 
deflned microscopic model of a Bose-condensed gas, and as a result, the relaxation times 
responsible for the dissipation can be calculated explicitly. 



We briefly reported on some of our main results in Ref. [|l^] which discusses the ZGN' 
equations for the special case of a uniform Bose gas. Here we give a more complete derivation 
of the ZGN' equations and a more detailed discussion of their properties. Our work may 
be viewed as an extension to the case of a trapped Bose-condensed gas of the pioneering 
work of Kirkpatrick and Dorfman |jl3| (referred to in the following as KD). However, these 
authors concentrated on the derivation of the Landau two-fluid equations with the inclusion 
of dissipation, and did not address the interesting possibility of the condensate and non- 
condensate components being out of diffusive local equilibrium {fic i^)- In this context, 
we should also mention Bogoliubov's rigorous fleld-theoretic derivation of the Landau two- 
fluid equations fl^ , which also implicitly assumed that the superfluid and normal fluid were 
always in local thermodynamic equilibrium. 

In the last few years, several papers dealing with non-equilibrium properties of trapped 
Bose gases have appeared. We mention speciflcally the work of Gardiner et al. [|T5|, Proukakis 
et al. Il^], Stoof and Walser et al. [Q. Each of these papers presents a formal derivation 
of a quantum kinetic equation using quite different theoretical methods. The equations of 
motion we derive for the condensate and non-condensate are closest in form to those of 



Stoof and Walser et al. [|18|, although there are also many formal similarities with the 
work of Proukakis et al. . However, in contrast to this other work, we focus on deriving 
hydrodynamic equations, and on establishing the criteria under which such equations are 
valid. In particular, we identify the crucial role that collisions between condensate and 
non-condensate atoms have in estalishing a state of local equilibrium. The paper of Stoof 



Tl\ is also useful for its thorough and insightful discussion of the dynamics of condensate 



formation and growth, a topic which Gardiner et al. [|I^ study in considerable detail. In 
this regard, we note that our kinetic equation can also be used to investigate condensate 



growth and applications to this problem will be presented elsewhere |19|] . 

A brief summary of our paper follows. In Section II, we derive the equation of motion for 
the condensate, and in Section HI, we obtain the quantum kinetic equation describing the 
non-condensate. These equations are applicable to the coupled dynamics of the condensate 
and non-condensate at flnite temperatures, and have been derived taking careful account 
of correlation functions involving three fleld operators, which describe the transfer of atoms 
between the two components. Taking moments of our kinetic equation, we derive in Section 
IV the hydrodynamic equations for the non-condensate, based on the assumption of local 
equilibrium within the excited-atom component. In Section V, we then analyze the full 
set of linearized ZGN' equations and show that a relaxation time naturally arises which 
characterizes the equilibration of the condensate and non-condensate chemical potentials. 
We show (in the simple case of a uniform Bose gas) that the ZGN equations give a good 
description in the limit cur^ ^ 1, while the Landau two- fluid hydrodynamic theory emerges 
in the opposite limit, ojt^ ^ 1. In Section VI, we demonstrate that the generalized Gross- 
Pitaevskii equation for the condensate and the kinetic equation for the non-condensate are 
consistent with the generalized Kohn theorem for the center of mass motion in a harmonic 
conflning potential. In Section VII, we reformulate the ZGN hydrodynamic equations in 
terms of a variational functional which is then used to calculate the frequencies (and as- 
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sociated condensate and non-condensate velocities) of the lowest-lying monopole (/ = 0), 
dipole (/ = 1) and quadrupole (/ = 2) collective modes. The results for the out-of-phase 
dipole mode were first reported in ZGN |0. Finally, in Section VIII, we assess what we have 
accomplished and briefly discuss some further extensions of the theory. 



II. HYDRODYNAMIC EQUATIONS FOR THE CONDENSATE ATOMS 

We start with the usual Heisenberg equation of motion for the quantum field operator, 

,n^}t^ = [^(r, t),H]= (^-^ + UUr)^ ^(r, t) + #t(r, t)^(r, t)^(r, t) , (1) 

where Uext{T^) is the confining potential. We have here assumed that the interaction potential 
can be represented as a short-ranged pseudopotential of strength g = Airah'^ /m, where a is 
the s-wave scattering length. The equation for the condensate wave function is obtained 
by taking an average of (|l|) with respect to a broken symmetry nonequilibrium ensemble in 
which the field operator takes a non-zero expectation value (denoted by angular brackets), 

$(r,t) = (V^(r,t)). (2) 

Introducing the usual definition of the non-condensate field operator ^/'(r,t) according to 

^(r,t) = $(r,t)+V^(r,t), (3) 

with (?/'(r,t)) = 0, the expectation value of (|l|) thus yields 

' = ( - ^ + Ue,t{r) + gn,(v, t) + 2gfi{v, t) \ <l>(r, t) 

+ gmiv, t)$*(r, t) + g{^\v, t)ij{v, t)^{r, t)) , (4) 

where nc(r,t) = |$(r, t)p is the non-equilibrium density of atoms in the condensate and 
h{r,t) = {ilj''{r,t)ilj{r,t)) is the non-equilibrium non-condensate density. Eq. @ also in- 
volves the off-diagonal non-condensate density m(r, t) = {ilj{r,t)ilj{r,t)) and the three-field 
correlation function (ip'^tptp) , both of which in principle have non-zero expectation values 
because of the assumed Bose broken symmetry. We note that is formally exact within 
the pseudopotential approximation. A more detailed discussion starting from the bare pair 
interaction is given in Ref. . 



It is useful at this point to discuss some standard approximations to (^. It reduces to 
the well-known Gross-Pitaevskii (GP) equation if all the atoms are in the condensate 
(i.e., n = 0) and the anomalous correlations (m and {ip'^ipip)) are ignored. This is a very 
good approximation for T ^ Tbec] at T = 0, the non-condensate fraction in trapped 
atomic gases is estimated to be less than 1% pO|,3. As a first step in dealing with finite 
temperatures, Hutchinson et al. (referred to in the following as HZG) calculated the 
excitations in the HFP approximation |2T| . The thermal occupation of these quasiparticle 
excitations provides an expression for the equilibrium non-condensate density no(r), while 
the condensate density nco(r) = |$o(r)P is obtained from the solution of a generalized GP 
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equation. The self-consistent calculation of these two densities determines the equilibrium 
properties of the trapped gas within the Popov approximation and the excitation energies 
are themselves identified with the collective excitations of the condensate. In the context 
of Eq. the HZG approximation consists of neglecting m and {ip'^ipip) and solving the 
linearized version of the resulting equation with the additional approximation n(r, t) ^ no{r). 
In physical terms, this latter approximation corresponds to the condensate atoms moving in 
the static Hartree-Fock field due to the non-condensate atoms. That is, only the dynamics 
of the Hartree mean field of the condensate atoms is treated self-consistently, while the 
collective dynamics of the non-condensate atoms is completely ignored. 

Returning to the exact equation for the macroscopic wavefunction in (^), it is useful to 
introduce phase and amplitude variables, $(r, t) = |$(r, t)|e*^*- After a little algebra, one 
finds that is equivalent to 

^ + V ■ Kv,) = ^Im[(<|.*)2m + <I>*{^^M (5a) 

= - ^{V0f = -(/., + ^mv'^). (5b) 

Here the new local variables are the condensate density 

n,{r,t) = mr,t)f, (6) 



and the condensate velocity 



v,(r,t) = Ave(r,t), (7) 

m 



while the local condensate chemical potential is defined by 



h^Vyn,iT,t) 

/ic(r, t) = . + Uext{r) + gricir, t) + 2gn{r,t) 



+ -^Re[(<l>*fm + <l>*(i!^iP)]. (8) 
nnc 

The condensate velocity in our theory is clearly identifiable with the superfiuid velocity in 
the usual discussions of superfiuid hydrodynamics. It will be seen later when we discuss the 
kinetic equation for the excited atoms that, in our model, the condensate density Uc can be 
identified with the superfiuid density as defined in the Landau theory 

One may give a physical interpretation of the generalized condensate equation in (|5b| ). 
In static thermal equilibrium, we have = and all variables are time-independent. In 
this case, the RHS of ( ^a]) must vanish and (|5b|) has the solution 6{t) = —ficot/fi- Thus it 



is clear that in thermal equilibrium, a condensate atom has an energy Sco equal to jj^co, the 
equilibrium chemical potential. We now suppose that we are in a regime where the variables 
He and Vc are slowly-varying in space and time. In this situation, we have 

96 

e{r, t) ~ ^(ro, to) + ^(t - to) + ■ (r - ro) + . . . 

1 1 
= 9{to, to) - ^^c(ro, to){t - to) + -mVc(ro, to) ■ (r - fq) + . . . , (9) 
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where 5c(ro,^o) = Aic(ro,^o) + \^v^i''^o,to)- With this phase variation of the condensate 
wavefunction, it is natural to interpret Ec as the local energy of a condensate atom and 
Vc as its local velocity. This interpretation is consistent with the energy and momentum 
conservation laws we derive later when we consider collisions between condensate and non- 
condensate atoms (see, in particular, the discussion at the end of Section ^ and in Appendix 
A). 

The condensate equations summarized by (|^-^) are formally exact, but are not a closed set 
of equations since they depend on the additional non-condensate variables n, m and {ip'^ipip) . 
As in Ref. [0, we restrict ourselves in the present paper to the Popov approximation |^ 



which corresponds to neglecting the anomalous pair correlation rh in (^) and (g). Some 
justification for this approximation is given in Appendix A. However a crucial difference 
from ZGN is that we retain the three-field correlation function (ip'^ipip) in (^a|) and which 
we evaluate in an approximation consistent with setting m to zero. We should remark that 
the importance of {ip'^ipip) when dealing with collisions involving condensate atoms has also 
been emphasized in Ref. [p!6|] . 



III. DERIVATION OF A QUANTUM BOLTZMANN EQUATION 

In this section, we turn to the dynamics of the non- condensate. The physical properties 
of interest are in principle defined by the following equation of motion obtained from (|l]) 
and d), 

QJ, ( ^2y2 \ 

ifi-^ = I - + Ue^t + 2gnU- 2ghiP + g^'^^P^ 

+ g^^iijjip - m) + 2g^{i)'^^i> - n) + gi'ip^^} - {i>''i>i>)) , (10) 

where n = ric + n is the total density. This equation preserves = as a function of time. 
It will allow us to derive a kinetic equation for the excited atoms. Following Kirkpatrick 



and Dorfman |]T3[, it is convenient to define the time evolution of ip by 

iP{r,t) = S^{t,to)i'{r,to)S{t,to), (11) 
where the unitary operator S{t, to) evolves according to the equation of motion 

th^^^ = H^,it)Sit,to), (12) 

with S(to,tQ) = 1. Here, to is the time at which the initial nonequilibrium density matrix 
p{to) is specified. The effective Hamiltonian in ([T^) is given by 



H'it) = H[{t) + H',{t) + H',{t) + H'S) , (13) 



where the various contributions are defined as 
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Ho{t)= J dr^U~^ + U{r,t)U 
H[{t)= Jdr (Li(r,t)^t + ^*(r,t)^) 

H'^{t) = ly dr^^^^ ~ d'cn{v,t)'4)^ 
f/(r, t) = Ue.t{r) + 2gn{r, t) = U,.,t{j) + 2g[fi{j, t) + n,(r, t)] 
Li(r, t) = -g (2n{r, t)<l>{r, t) + m(r, t)$*(r, t) + (^^V^^)) . (14) 



It is understood that in all these expressions the arguments of the field operators are (r, to); 
the dependence of the various terms of the Hamiltonian on t arises through the mean-field 
expectation values. It can be shown that Eqs. (|1T])-(|T^, together with the equal time 



commutator [?/^(r, tg), "^/'^J''' ^o)] = ^(r ~ r'), reproduce the equation of motion for ip in (|T0|), 
as well as that for ip'^ . In writing (|TBp we consider H'{t) to be a perturbation to the zeroth- 
order Hamiltonian Ho{t). Noting that U{r,t) is the total Hartree-Fock (HF) mean field, 
Ho{t) defines excitations of the system at the level of the time-dependent HF approximation. 
Other choices of the zeroth-order Hamiltonian are also possible. For example, Ho{t) could 
be combined with H2{t) to define a dynamic Bogoliubov Hamiltonian which would then lead 
to local Bogoliubov quasiparticle excitations. This possible extension which would be 
appropriate when considering very low temperatures, will not be considered in this paper. 

We now consider an arbitrary operator 0{t) which is made up of some combination of 
non-condensate field operators ip{r,t) and '?/'^(r, t), for example, the local density operator 
h{r,t) = ip''' {r,t)ip{r,t). By making use of (0), the expectation value of 0(t) with respect 
to the initial density matrix p(to) can be expressed as 

(6(t)) = (O), = Trp(to)0(t) 

= Trp(t,to)0(to), (15) 

where pit, to) = S(t,tQ)p(to)S'^ (tjto) satisfies the following equation 

'^^^^ = [HeE{t)rp{t,to)]. (16) 

We refer to Appendix A for further discussion of this equation of motion. 

Our ultimate objective is to obtain a quantum kinetic equation for the non-condensate 
atoms. For this purpose, we define the Wigner operator 

/(p,r,to) = / rfr'e^P-'/?i^t(r + ^,to)^(r-^,to). (17) 
Its expectation value according to Eq. ([T5|) then yields the Wigner distribution function 

/(p,r,t) = Trp(t,to)/(p,r,to). (18) 
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Knowledge of this function allows one to calculate various nonequilibrium expectation values, 
such as the non-condensate density 

Using (|16D , the equation of motion for / is 

^^^^^ = ^Trp(t,to)[/(p,r,to),^eff(t)] 

= ^Trp(t, to) [/(p, r, to), ^o(t)] + ^Trp(t, to) [/(p, r, to), H\t)] . (20) 

The first term on the right hand side of (^) defines the free-streaming operator in the kinetic 
equation. With the assumption that U (r, t) varies slowly in space, we then have 

+ £ ■ V/(p, r, t)-VU- Vp/(p, r, t) = J-Trp(t, to) [/(p, r, to), H'{t)] . (21) 

The right hand side of this equation clearly represents the effect of collisions between the 
atoms. The reduction of this term to the form of a binary collision integral is a lengthy 
exercise (some details are given in Appendix A). However, the final result has a physically 
transparent form. The collision integral is the sum of two contributions: 



dl 

dt 



= Cu[f]+C22[f], (22) 

coll 



where 



C22[f] = 7^%^ / ^P2 / ^P3 J dpiS{p + P2 - P3 - P4) 



(27r)^ 

x5{e, + Bp, - Sp, - 6,,) [(1 + /)(1 + /2)/3/4 - //2(1 + /3)(1 + U)] , (23a) 



Ci2[f] = ^c^^Y^i J ^Pi / ^P2 / dp36{mvc + Pi - P2 - Pa) 



x6{ec + Ep, - Ep, - Ep,)[6{p - Pi) - 6{p - P2) - 6{p - ps)] 

X[(l + /l)/2/3 - /l(l + /2)(1 + /3)], (23b) 



with / = /(p, r, t), fi = f{pi,r,t). We also note that the variables Uc, Vc, Ec and ip are 
all functions of r and t. As a result of Bose statistics, the creation of an atom in a state i 
is associated with the statistical factor (1 + fi). Clearly C22 describes two-body collisions 
between excited atoms (2 atoms ^ 2 atoms) while C12 describes collisions which involve 
one condensate atom (1 atoms ^ 2 atoms). The momentum and energy delta functions in 
the C12 collision term take into account the fact that a condensate atom locally has energy 
Ec{r,t) defined in (^), and momentum mvc. On the other hand, the excited atoms locally 
have the HF energy 
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2 

ip{r,t) = ^ + U{r,t). (24) 

We note that the HF mean-field defined in involves the condensate density nc{r,t) = 
|$(r,i:)p, which must be determined self-consistently using the condensate equation of mo- 
tion in (g). 

To complete our derivation of a closed set of equations for both the condensate and 
non-condensate, we must finally deal with the off-diagonal HF self-energy m and the term 
(ip'^ipip) which appear in and (j^). We show in Appendix A that m can be consistently 
neglected, and that 

(iP^iP) = J dpi J dp2 J dp36{mY^ + Pi - P2 - Ps) 

x5{e, + e,, - e,, - e,,)\h{l + h){l + h) - (1 + /O/a/s]- (25) 

We do not display the real part of {ip'^ipip) which contributes to the local condensate chemical 
potential defined in (§). This contribution is of order and will be neglected since our 
approximations do not adequately treat all terms of this order (see Appendix A for further 
discussion). However, the imaginary part of {ip^ipilj) shown in (|25|) contributes an important 
source term, — ri2[/], to the right hand side of the continuity equation in (^). Comparing 
(^) with the expression for Ci2[/] in ( |23b|) , this source term is given by 

- ri2[/] = |lm[$*(^W)] = -/ (^C'i2[/(p,r,t)] . (26) 

The latter form in (^) emphasizes the fact that the Cu collisions can locally change the 
relative number of condensate and non-condensate atoms and that separate number con- 
servation laws no longer hold for each of these components. As we shall see, the close 
connection between {ip^ipip) in ( pS]) and the C12 collision integral in ( P3b|) is essential for the 



conservation of the total number of particles. This indicates that the approximations leading 
to our final equations for the condensate and non-condensate are internally consistent. 

To summarize the results of this section, as well as that of Section we write out the 
final form of the equations we have derived for the dynamics of a trapped Bose-condensed 
gas at finite temperatures. For the condensate, we have the equations 

^^^ + V-(n,v,) = -ri2[/] (27a) 



dt 
'd 
dt 

where the condensate chemical potential is given by 



m(-^ + v,- v)ve = -V/i,, (27b) 



h^V\nJr,t) 

/i,(r, t) = y + Ue^t{r) + gn,{r, t) + 2gfi{v,t) . (28) 

2mWnc(r, t) 



As we have discussed in detail, these equations are approximations to the exact equations 
given in (Ef) and (H). For the non-condensate we have the kinetic equation 
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^^^^r^ + - ■ V/(p, r, t)-VU- Vp/(p, r, t) = C^,[f] + C^^lf] , (29) 
at m 



with the two colhsion terms given exphcitly by the expressions in (|23D. These equations 
are coupled through the HF mean fields and must therefore be solved self-consistently. An 



equation equivalent to (^) was derived by Proukakis et al. fl^, while our kinetic equation 
in (p9|) is implicitly contained in the work of Stoof |1T7| . 

It is apparent that the equations in (^) are already in the form of "hydrodynamic" 
equations for the condensate variables. In Section 0, we use (^) to derive analogous hy- 
drodynamic equations for the non-condensate atoms. Although these equations have been 
derived for a trapped gas at finite temperatures, we note that they go over smoothly to the 
correct T = theory, since as /(p, r, t) and hence n become negligible, the condensate equa- 
tions ( P7D and ( pSD reduce to the standard Gross-Pitaevskii dynamics of a pure condensate 
at T = [|l|,0. We therefore expect our equations to provide a reasonable description of 
the dynamics over the full range of temperatures from the strongly degenerate limit well 
below Tbec through to the classical high-temperature limit. 

Before leaving this Section, we compare these results with some earlier theoretical stud- 
ies. For temperatures above the Bose-Einstein transition (T > Tbec) where C12 vanishes, 
the kinetic equation (|29|) reduces to the well-known Uehling-Uhlenbeck equation pSf , as de- 



rived in Ch. 6 of Kadanoff and Baym for example, and studied numerically by several 
authors [^. Below Tbec, (|29|) is analogous to the kinetic equation derived by KD for 



a homogeneous weakly-interacting Bose-condensed gas at finite temperatures. The main 
difference is that these authors choose to work with excitations which are defined within the 
local rest frame of the superfluid (condensate), and thus the excitation energies are measured 
relative to that of the condensate atoms. In contrast, our kinetic equation involves atoms 
moving in a self-consistent HF field. It is perhaps useful to make this comparison more 
precise. Our Ci2[/] collision integral involves the energy- momentum conservation factors 

5{ec + 61-82- es)S{mVc + Pi - P2 - Ps) 

where ii = pf/2m + 2g{n + ric) (assuming Uext = for simplicity) and Ec = fic + \'mv1 is the 
local energy of an atom in the condensate. On the other hand, in the work of KD, Ci2[/] 
contains the energy- momentum conservation factors 

5{E, -E2- Es)S{p[ p's) , 

where Ei is the quasiparticle (Bogoliubov) excitation energy and p^ = Pi — is the 
quasiparticle momentum in the local rest frame. In the "high-temperature limit" (large 
momenta), these excitation energies reduce to E^ ~ pf/2m + gric = ii — Hc — v ■ pi + |mw^ 
where Hc = gi^c + 2ri) is the local condensate chemical potential as given by (|2^) . One 
can easily verify that the two energy-momentum conservation factors defined above are 
completely equivalent in this high-temperature limit, and thus the two seemingly different 
formulations are essentially the same. 

It should be noted that the high-temperature limit in the case of a trapped Bose gas 
is not as restrictive as it might seem. The work reported in Ref. indicates that the 
differences arising from using the HF vs. Bogoliubov excitation spectrum in calculating 
thermodynamic properties are very small down to quite low temperatures, in contrast to 
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the situation for a uniform Bose gas. A more detailed explanation for the dominance of 
single-particle excitations in determining the statistical behaviour of trapped Bose gases can 
be found in Ref. 



IV. NON-CONDENSATE HYDRODYNAMIC EQUATIONS 



A. General moment equations 



Following the standard procedure in kinetic theory, we proceed to take moments of the 
kinetic equation in (pQ]). It follows immediately from their definitions that 



dp PC22 = 0, J dp {p- mWc)C_ 
dp epC22 = 0, j dp {ip - ec)C 



12 = 
12 = 0. 



(30) 



This shows that both the Cu and C22 collision processes conserve energy and momentum. 
In the case of 1 # 2 collisions, the momentum of the non-condensate atoms is conserved in 
the local rest frame of the condensate; in addition, the energy conservation condition takes 
into account explicitly the different mean-field energies of condensate and non-condensate 
atoms. In addition to the exact relations in fBOl), we have 



j dp C22 



0, 



(31) 



which implies that the number of excited atoms is conserved in C22 collisions. In contrast, 
as already mentioned after (^), C12 does not conserve the number of excited atoms and 
thus we find 



ri2[/] 



dp 



C12 ^ 0. 



(32) 



Taking (^)-(|32|) into account and, after a certain amount of rearrangement, the exact 
moment equations derived from ( p9D can be written in the form {fi and u denote Cartesian 
components): 



^ + V-(nv„) = r,2[/] 



(33a) 



~ d \ dP,, ^ dU 

mn \ — + \n ■ \ Vn^i = TT^ U — 

\Ot J OXy OXfj, 



- n- m{vn^, - Vcf,)Ti2[f] 



(33b) 



^ + V.(ev„) 



-V ■ Q - D^yP^y 



1 , 



ri2[/]. 



(33c) 



The non-condensate density is defined in (0) while the non-condensate local velocity is 
defined by 
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,Xr,iK(r,i)-/p^|p£/(p,r,i). (34) 



In addition, we have 



(27r^)'^ \m 

^ (2^2^^"^ ~ V(p, r, t) . (35c) 



The symmetric rate-of-strain tensor appearing in (33c) is defined by 




It is clear from the addition of (|27a|) and (|33a|) that the total number of atoms is conserved, 
as it must be. 



The set of equations in (p3D -(p5|) are formally exact but they are obviously not closed 
since the number of local variables (n, ric, v^,, as well as the stress tensor P^,^, kinetic 
energy density e, heat current Q, and ri2) exceeds the total number (nine) of coupled scalar 
equations. To proceed we must specify the conditions under which the dynamics of the 
trapped gas is to be determined. The range of possibilities include: (i) the collisionless regime 
in which C12 = C22 = 0. Here, the collisionless kinetic equation for the non-condensate 
must be solved explicitly in order to determine the various physical quantities of interest; 
(ii) an intermediate regime in which collisions between excited atoms are sufficiently rapid 
to establish a state of local equilibrium within the non-condensate component, but in which 
collisions between the condensate and non-condensate, as described be Cu, can be neglected. 
We shall refer to this as a state of partial local equilibrium which might be expected to 
arise near Tbec when the density of condensate atoms is small (note the proportionality of 
C12 to ric); and (iii) a regime in which the collision rates are so high that a condition of 
complete local equilibrium is estabhshed. As will be discussed in detail later, regime (ni) is 
the one conventionally dealt with using Landau two-fiuid hydrodynamics in uniform Bose 
superfiuids 0. Of course, the most general situation will not fall neatly into any of the 
categories listed so far and will require a detailed solution of the kinetic equation in (^9[) 
together with the quantum hydrodynamic equations for the condensate. 

We shall not deal with the collisionless regime (i) in the present paper, but note that 
some work along these lines has recently appeared ||2^. Instead, we shall focus on the other 
two regimes indicated above, in which the C12 and C22 collisions play a crucial role. Of 
particular interest is the transition from the partial to complete local equilibrium conditions. 
An analysis of this regime will clarify the way in which the Landau theory emerges within 
the context of our explicit microscopic model. This turns out to be a rather subtle problem, 
which was not fully appreciated in our earlier work P,p8|, nor addressed in the analogous 
work of KD 
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B. Partial local equilibrium 



We now consider the regime in which Eqs. reduce to a closed set of equa- 

tions for the non-condensate variables. When the collision rate among the excited atoms is 
high, the collision integral C22 drives the distribution function /(p,r,t) towards the "local 
equilibrium" Bose distribution 

= e/3[^(p-mv„)2+C/-A] _ 1 ' 

where now the temperature parameter /3, local fluid velocity v„, chemical potential /i and 
mean field U are all functions of r and t. It is important to appreciate that the local chemical 
potential fl which appears here is distinct from the local chemical potential /ic defined in 
(p8|) for the condensate. How these two chemical potentials acquire a common value in the 
limit of complete local equilibrium requires a careful analysis of the effect of Cu (see Section 
0). For the moment, we shall simply note the consequences that ( ^7|) has for the collision 
integrals. 

First, one may immediately verify that the local equilibrium form of / in (l37| ) guarantees 
that 

C22[f] = 0. (38) 

Indeed, it is this condition that defines the "local equilibrium" solution. The result is 
independent of the value of fl and only makes use of the energy and momentum conservation 
factors in the definition of C22 in ( P3a| ), and the key identity for the Bose distribution 

l + f{x) = ^^ = -f{-x). (39) 
e — i 

In contrast, one finds from ( |23b|) that Ci2[/] 7^ 0, namely 

X J dpi J dp2 J dps5{mvc + pi - P2 - P3)5iec + ei - £2 - £3) 

X [5(P - Pi) - Sip - P2) - 5{p - P3)](l + /l)/2/3. (40) 

A factor similar to the first square bracket in (^) appears in the growth equation derived 
by Gardiner et al. [l^. In their work, our time- and space-dependent jl is replaced by a 
constant chemical potential characterizing a non-condensate particle reservoir. On the other 
hand, they use a kinetic description for low-lying non-condensate states which we describe 
in terms of a local equilibrium distribution with local chemical potential jl. Thus, although 
different in detail, there is considerable overlap in physical content. The important point 
is that Ci2[/] is in general non- vanishing, reflecting the fact that excited atoms in the non- 
condensate need not be in local diffusive equilibrium with the condensate atoms [jl 7^ /Xc)- 

Assuming that the collision rate among excited atoms is sufficiently high, the local equi- 
librium distribution in (|37D provides the appropriate zeroth-order solution to the kinetic 
equation in (^). This is the same reasoning used in classical kinetic theory ||2^. Of course. 
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even in the absence of the C12 colhsion term, / is not an exact solution of (|29D . A systematic 
procedure based on the Chapman-Enskog method gives rise to corrections to / which are 
necessary in order to describe thermal conduction and viscosity. A detailed discussion of 
these damping processes in a trapped Bose gas above Tbec is given in Refs. 13 and ^ 



Such calculations can be extended to below Tbec and have been carried out by KD [|T^ for 
a uniform Bose-condensed gas. However, we ignore these corrections in the present analysis 
and proceed with the assumption that / ~ / is a good approximation (see Section |V111| for 
further remarks.) 

Using ( P7| ) to evaluate the quantities in (P5D, we find that the heat current Q(r,t) = 0, 
and that 

P^(^. t) = S^P(r.t)^i,^lj^^f(p,r,t)l^. (41) 

In addition, we note that the kinetic energy density ( |35c| ) is given by 



e(r,t) = ^P(r,t), (42) 

which is the relation found for a uniform ideal gas. Since the momentum dependence of the 
integrand in (^Tj) is the same as for a static equilibrium Bose distribution, the integrations 
can be carried out explicitly to give the expression, 

Hr,t) = ^^g,/2{z) , (43) 

where the local thermodynamic variables (3 and z are again functions of r and t. The 
Bose- Einstein functions are defined as gn{z) = J^'iZi /l^, the local equilibrium fugacity is 

^(r,t)=e'^('-'*)['^('-'*)-^(^'*)l, (44) 

and the local thermal de Broglie wavelength is 



Finally, we note that the non-condensate density associated with ( pTf ) is given by 



To summarize, using the local equilibrium approximation / ~ /, the hydrodynamic 
equations in (^) simplify to 

— + V-(nv„) = ri2[/], (47a) 
mh + v„, • V^j v„ = - VP - nVU - m(v„ - v,)ri2[/] , (47b) 
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dP 
'dt 



+ V • (Pv„) 



-m V. 



n 



ri2[/], 



(47c) 



where 

ri2[/] 



(27r/i)3 



Ci2[~f] 



-/3(/i-/^c-5m{v„-Vc)^) 



(ipi / (ip2 / (ip, 



x5(mVc + Pl - P2 - P3)<5(£^c + ^1 - ^2 - e'3)(l + /l)/2/3 



(48) 



It is clear that the ri2[/] terms in (^7]) will play a crucial role in the solution of these 
hydrodynamic equations. In particular, these terms describe the transfer of atoms between 
the condensate and non-condensate and are responsible for bringing the two components 
into complete local equilibrium. 



In this paper, we concentrate on solving (|47|) for small amplitude oscillations around 
static (i.e., absolute) equilibrium. For the non-condensate atoms, absolute equilibrium is 
specified by a Bose distribution at a uniform temperature Tq and chemical potential fiQ (all 
static equilibrium quantities are indicated by a subscript 0). The equilibrium density and 
pressure are then given by the expressions fiQir) = g^/2{zo)/ and -Po(r) = 5'5/2(2;o)//3o^o; 
respectively. These are of the same form as in Eqs. (^61) and (^31) , but with all thermodynamic 
quantities taking on their equilibrium values. In particular, the equilibrium fugacity is 
zo{r) = eMi^o-Uoir)) ^^^^ Uo{r) = Uextij) +25^[no(r) +nco(r)]. It is clear that the equilibrium 
condensate density ncQ{r) is also required to complete the definition of the non-condensate 
equilibrium variables. This we obtain from the equilibrium condensate wavefunction which 
is a solution of the static equation 



(49) 



This is the HFP equation for the condensate used previously PJ26[], but with the important 
difference that no is determined in the HP approximation as described above. 
The eigenvalue of (HO) can be written as 



^'VV^co(r 



2mJnco(r) 



(50) 



and is identified with the static equilibrium condensate chemical potential. To see the 
connection between /iq and fico, we refer to the expression for Ci2[/] in Eq. (^Ol). With 
Vn = Vc = 0, we note that the factor in the square brackets vanishes if fiQ = fico, which 
is precisely the expected condition for diffusive equilibrium between the condensate and 
non-condensate. Thus the constraint /io = /ico must be imposed in order to obtain the self- 
consistent solution for no{r) and nco(r). Further discussion regarding the static equilibrium 
properties of a trapped Bose gas is given in Appendix B. 
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C. Linearized hydrodynamic equations 



To study the collective modes, we linearize the system of equations in ( P7| ) around the 
static equihbrium solutions given above to obtain the non-condensate hydrodynamic equa- 
tions 



mrio- 



d6n 

~df 
d6vn 
dt 
dSP 



-V ■ (no(5v„) + 6Ti2, 

-VSP - 5hVUo - 2ghoV{Sh + 5nc), 

5 ~ 2 ~ 2 

--V • (Po5v„) + -5v„ ■ VPo + -(/ico - Uo)5Tu 



(51a) 
(51b) 

(51c) 



We observe that 6T12 does not appear explicitly in ( |51b| ) as it contributes a second order 
correction. These non-condensate equations must be solved in conjunction with the following 
linearized equations for the condensate [see {^7\)] 



m- 



dSric 
dt 



where 



5/ic(r,t) 



1 



2<l>nW 



-6nc{r,t) 



+ gdricir, t) + 2gSn{r, t) . 



(52a) 
(52b) 

(53) 



The operator i2o(r) appearing in this equation is defined as 



Mr) 



2m 



+ Uext{r) + 2gno{r) + gnco{r) - jj^o 



(54) 



and is essentially the HFP Hamiltonian in Eq. (^9|). We must finally obtain an expression 
for 5ri2 in (BTf) and (E^) which is the linearized form of the expression in (HH). We have 



Pon. 



cO 



ri2 



diff 



(55) 



where we have introduced the difference between the local chemical potentials of the con- 
densate and non-condensate, 



/^diff = — l^c, 

and the "equilibrium" C12 collision rate 



(56) 



J dp, I dp2 I rfp35(pi - P2 - P3)5(/ico + e, - 62 - e,)(l + /D/aVs • (57) 



ri2 {2Tifh' 

Here, ff is the static equilibrium Bose distribution function. 
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D. ZGN hydrodynamics 



In the remaining part of this Section, we suppose that on the time scale of the collective 
modes of interest, collisions between condensate and non-condensate atoms can be ignored. 
This amounts to setting Ci2[/] (and hence ri2[/]) to zero, and gives the hydrodynamic 
equations derived in ZGN |Q. The linearized version of these equations is obtained by 
setting 6T 12 to zero in (^) and (|52|). They are a closed set of equations for the variables 
Suc, Sn, Vc, v„ and 6P and may be solved to obtain the self-sustained oscillations of the 
trapped gas at finite temperatures. A powerful variational method for solving equations of 
this kind is presented in Section [Vll] , where we consider the explicit solution of the ZGN 



equations for a trapped gas. 

An approximation commonly used in the recent literature on trapped Bose gases involves 
the neglect of the non-local "quantum pressure" terms, i.e., the first terms on the RHS of 
(0) and (|53|) . This so-called Thomas- Fermi (TF) approximation |2^j3^ has been found 



to give quite good results when the number of atoms in the trap is large. Within the TF 
approximation, the equilibrium condensate chemical potential in (pOD simplifies to 



/ico - f^ex't(r) + gricoir) + 2gno{r) . (58) 



At T = 0, where noi^r) = 0, ( pq ) gives a condensate density which goes smoothly to zero 
at points where t/ea:t(r) = /^co- However at finite temperatures, the combined solution of 
(0) and the semiclassical expression for the non- condensate density (no = (73/2 (2^0 )/Ao) 
gives rise to unphysical discontinuous behaviour of the densities at the sharp boundary of 
the condensate |Q. This behaviour is especially problematic when attempting to obtain 
solutions of Eqs. (|51|)-(|53D for the collective modes, since elaborate boundary conditions 
must be formulated in order to join the fluctuating variables at the condensate boundary. 
These problems can be avoided by retaining the quantum pressure terms which ensure that 
the densities are smooth functions of position in all regions of the trapped Bose-condensed 
gas. Although these gradient terms might appear to complicate the analysis considerably, 
it will be seen that they can be handled easily using the variational method of solution 
developed in Section |V11| . 

For the purpose of comparison with results obtained in Section 0, it is useful to consider 
the ZGN equations in the limit of a uniform Bose gas. In this case, the ZGN equations 
yield first and second sound modes p8[. Eliminating the variables 6nc, Sn and 6P from 



Eqs. ([5T|)-(|5^ (with 5T 12 = 0), one obtains two coupled equations pS 



^~d^ = gncoV{V ■ 5ve) + 2ghoV{V ■ 5v„) (59a) 

m^-^ = 2gn,oV{V ■ 5v,) + f + 2gnoj V(V ■ 5v„) . (59b) 

Introducing velocity potentials according to 5vc = V0c and 5v„ = V0„ and assuming a 
plane wave solution of the form 0c,n(i',^) = <Pc,n^^^^'^~'^^K one obtains 

mu^cpc = gn^ok'^cpc + 2ghok^(l)n (60a) 

mwVn = 2gncok'^(j)c + ( |^ + 2gho] k'^cp^ ■ (60b) 

V 3^0 / 
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The nontrivial solution of these equations yields the normal mode frequencies of the uniform 
gas. These modes are first and second sound waves with the dispersion relation uo = uk, 
where the sound velocities u are the solution of 



u -u [-^ + + + [-— =0. (61) 

Writing this equation as u'^—Au'^+B = 0, we see that the coefficients A and B are completely 
determined by equilibrium properties of the gas. These coefficients will be compared in 
Section with the corresponding coefficients obtained on the basis of the Landau two-fluid 
equations. As we show there, the A and B coefficients are slightly different in the Landau 
regime and reflect the different conditions underlying the propagation of flrst and second 
sound waves in the two cases. 

We note that the quartic equation in (|6TD yields four solutions while Eqs. (0) and (|52D 
are a total of flve equations and must therefore admit flve distinct longitudinal modes. The 
missing solution is, in fact, a zero frequency mode with mode amplitudes S\c = f^v^ = 0, 
6nc = —26n and SP = 2gnQSn. The physical interpretation of such a static mode is not 
immediately obvious. However, since we require some results from Section |V] in order to 
explain its meaning, we shall defer further discussion until then. 



V. COMPLETE LOCAL EQUILIBRIUM: LANDAU TWO-FLUID EQUATIONS 

In this Section, we show that within the TF approximation, the set of equations in 



and (]47|) can be used to derive linearized two-fluid equations which reduce in a certain limit to 
the standard Landau form 0, as generalized to include an external trapping potential. This 
derivation is of interest for several reasons. First, it shows that our approximate microscopic 
model is consistent with the Landau two-fluid equations when conditions of complete local 
equilibrium exist. Under these conditions, the Landau equations are valid for dense fluids 
such as superfluid ^He as well as for dilute Bose gases. These equations have been derived 



using very general arguments fTj], but this very generality tends to obscure the underlying 
physics which a microscopic model such as ours can reveal in a very clear manner. For 
example, our model not only leads to the Landau equations, but at the same time gives 
the static equilibrium thermodynamic functions required in the calculation of the normal 
modes. An even more signiflcant advantage of our microscopic model is that we can explore 
the transition from partial to complete local diffusive equilibrium. In this respect, our two- 
fluid hydrodynamics based on Eqs. ( P?] ) and (0) provides a more complete description of 
the possible dynamical behaviour in trapped Bose gases. 

In a previous paper we indicated how the linearized ZGN equations could be com- 
bined to give the Landau form of the two-fluid equations for a uniform gas. However, the 
arguments presented there were incomplete. Here, we present a more careful discussion 
starting form Eqs. (|27| ) and (^) and show how an explicit consideration of the source term 

[/] associated with the C12 collisions leads to generalized two- fluid equations that reduce 
to the Landau equations in the appropriate limit. In addition, we shall extend our previous 
analysis to the case of an inhomogeneous system. In doing so, however, we make use of the 
Thomas-Fermi approximation which effectively treats the system as locally homogeneous. 
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If this approximation is not made, one cannot obtain the usual Landau equations since 
the gradient term of the chemical potential in (|28|) contributes to the right hand side of 
( P7b| ). This nonlocal term is a quantum correction which is always ignored in microscopic 
derivations of the Landau equations for uniform systems [l^ . 



For completeness, we first display the linearized Landau two-fluid hydrodynamic equa- 
tions [|6|,|7|,|3^ , in the form they take in the presence of an external potential 



^ + V-5j = (62a) 



= -V6P - 6nVUe,t (62b) 



= -VSfi (62c) 

85 s 



dt 

where 



+ V-(so5v„) = 0, (62d) 



m(5j(r, t) = p,o(r)<^Vs(r, t) + p„o(r)(5v„(r, t), 

'm6n{r,t) = 6ps{r,t) + 6pn{r,t). (63) 

We have here distinguished the superfiuid variables ps and from the corresponding vari- 
ables for the condensate; a similar distinction applies to the normal fluid variables as opposed 
to the non-condensate variables. However, the condensate and non-condensate variables can 
be identified with the more usual two-fluid variables within the context of our model which 
treats the excitations in the HF approximation. It should be noted that this correspondence 
will no longer be valid if one were to go beyond this approximation. For example, in the HFP 
approximation, there is a finite depletion of the condensate even at zero temperature P] and 
the superfiuid density is no longer precisely equivalent to the condensate density. The other 
local variables appearing in the above equations are the pressure P, entropy density s and 
thermodynamic chemical potential p. As we shall show, each of these quantities is given by 
an explicit expression within our microscopic model. 

The linearized ZGN' continuity equations are given in ( |51a| ) and (|52a|) . Adding them. 



we obtain the continuity equation ( |62a| ). Combining our two velocity equations in ( 51b ) and 
( PB|) , we find 



— = -V5P - 2g{nco + no)V5h - g{n^Q + 2fio) V5nc 
ot 

- 2gV{nco + no)Sh - VUextSh . (64) 



One can show that (p^ is equivalent to (|62b|) if the total pressure is defined as 

P = P + ^gln"^ + 2nh - h^]. (65) 



21 



To prove this, we note that (|65D gives 

5P = 6P + 2gno6h + gijicQ + 2no)5nc 

and hence, 

V(5P = V5P + 2gnQV5h + g{ncQ + 2no) V5nc 
+ 2g'VnQ6h + gCVrico + 2Vno)<^^c • 



(66) 



(67) 



Making use of the TF result in (^) to rewrite the last term as —SricVUext, and inserting the 
resulting expression for V5P into (|62b|) , we obtain an equation identical to (|6^) . The results 
in (|6^-(|67D were given previously by ZGN for a uniform gas, but the present derivation 
is more general in that it is based on the ZGN' equations in the context of a nonuniform 
system. 

In order to make contact with the remaining Landau equations ( |62c|) and (|62d| ), we must 
introduce the appropriate thermodynamic variables corresponding to the chemical potential 
/i, local temperature T and the local entropy s. Within the HF approximation for the 
non-condensate, the total energy density is given by 



^+^9 (i^HHi^) + nUext 



;9 



n + 2nn 



+ nU, 



ext ) 



(68) 



where e is defined by (^) and (^21). For the local entropy we use the definition 

dp 



1 
T 



{27ihY 



;i + /)ln(l + /)-/ln/ 



Ip - n{fi - U) 



(69) 



This result corresponds to the entropy of a uniform ideal gas above Tbec with a fugacity 
z = e^^^~^\ In static equilibrium, (B^) reduces to 



5 ~ 

sqTq = -Pq + gucono 



(70) 



where we have noted that /to = /ico = Uext + 5'(2no + rico) and Uq = Uext + '^gno- 

From these expressions for P, e and s (all of which depend on r and t), we obtain the 
relation 



e + P — sT = ^c^c + fin = ficU + fidisn , (71) 

where fi^is is the chemical potential difference defined in (^61). As we shall see, this variable 
plays a crucial role in understanding the approach to complete local equilibrium. If fi^m were 
zero, then (^) would correspond precisely to the usual thermodynamic relation involving 
these variables, with playing the role of the equilibrium chemical potential. 

Having defined these various thermodynamic functions, we can now consider their local 
variations from absolute equilibrium (recall that the equilibrium properties are still functions 
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of position). We begin with the kinetic pressure given by (^31). Its variation leads to the 
equation 



6P = 6 



1 



= so5T + no{5jl-2g6n) (72) 

In obtaining this result, we have made use of ([70| ) and the Bose identity zg^{z) = gn-i{z). 
Combining (|66D and (|7^) , we find 

6P = so5T + no5/ic + no5 fidis ■ (73) 

The appearance of the last term is once again due to the non-condensate having a chemical 
potential different from that of the condensate. Similarly, from (|69D and (|7^), we find the 
variation of the entropy to be given by 

n6s = pP + gn,oSn. (74) 

This expression finally allows us to obtain an entropy conservation equation. Taking the 
time derivative of (|7^ , and making use of the non- condensate density and pressure equations 
( pla| ) and ( plcj ), we find 

d6s 3 85 P guco d5n 
"dt ~ 2To~dr ^ To dt 

= -V ■ (so5v„) . (75) 



We note that the source terms on the right hand sides of (|51a)) and ( plc| ) cancel out, so that 
a strict conservation law for the nonequilibrium local entropy defined in (^) is obtained. 

Finally, a comparison of ( |52b| ) and ( |62c[ ) shows that these two equations are equivalent 
if 6\s is identified with and 6fi with 6fic- With this identification, we arrive at a set of 
hydrodynamic equations which are precisely of the form of the Landau two-fluid equations. 
At first sight, it might therefore appear that our starting equations are in fact equivalent to 
the Landau equations, but this is not the case. The appearance of fi^is in (|7lD implies that 
He is not related to the other thermodynamic variables in the same way that the equilibrium 
chemical potential /x would be. As a result, 6fic and 6s are not related directly to the 
thermodynamic fluctuations 6T and 6P, as assumed in the usual derivations of the first and 



second sound velocities 34 



To see this more clearly, let us use (|7D|) to define the equilibrium local entropy function, 
Ts'^'^ = |P + gUcU, where the notation s'^'^ is used to distinguish this quantity from the 
local entropy given by (pP|). It is the variation of the equilibrium entropy s'^'^ which can be 
expressed in terms of the thermodynamic fluctuations 6T and 6P via equilibrium thermo- 
dynamic derivatives. By considering a quasistatic change in the thermodynamic state of the 
system, (^) leads to the relation 

5 ~ 

ToSs""^ = -5P + gricoSn + gho5nc - Sq^T, (76) 



and using ( [72D to eliminate 6T, we find 
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3 ~ 

TaSs""^ = -SP + gUcoSn + no5(Xdis ■ (77) 

This differs from the result for 5s given in {\f^. The last term in ([77| ) indicates that the 
variation of the equilibrium entropy is not simply given in the present situation by the 
equilibrium variation. In terms of ^s^*^, the entropy equation analogous to ( [75| ) takes the 
form 

= -V ■ (so5v„) + no—g^ . (78) 



The last term in ([78|) can be interpreted as the production of entropy associated with the 
equilibration of the condensate and non-condensate chemical potentials. Its appearance 
implies that the thermodynamic fluctuations 6P and 6T are coupled to fluctuations in 6^dis, 
which is directly related to the condensate- non-condensate collision terms appearing 
on the right hand side of Eqs. ( |51a| ), ( plc| ) and ( ^2a| ). To this point we have not made 
any assumptions about 5ri2, and as a result, the equations as they stand differ from the 
usual Landau two-fluid theory. To make this final connection we must obtain an equation 
of motion for 6fidis- 

Noting that ft — 2gn = fidis — gn-c (0) can be rewritten in the form 

dP = SoST + no{Sfidis ~ gSn^) ■ (79) 
Similarly, the variation of (|46|) yields the equation 

5h = (^ho + 7onco) ^ + — (5/idifr - gSuc) , (80) 
\/ / lo g 

where we have introduced the dimensionless quantity 

70 = ^gi/2izo) . (81) 



This quantity can be related to thermodynamic derivatives by noting that ( pOD implies the 
equilibrium variation at constant temperature, ShlT = {%/ g)6{fl — 2gn)\j' = — 7o(5n — 5fi) 1^. 
We thus conclude that 



'dn J 



T 



In a similar way, one can show that the coefficient of 5T in (80) can be expressed as 



^no + 7onco) ^ = (1 - ^o) • (83) 

Eqs. (^) and ( PDD show that the variations in 5T and 5/idifr can each be expressed 
in terms of 5?t,c, 5h and 5P. Eliminating from these equations and making use of 

Eqs. (|51a|) and ( |51c|) , we obtain the equation of motion for 5T, 
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1 d5T 2 g{no + llon,,) 

- V ■ dv„ - ~ " 6Tu ■ (84) 




^0 dt 3 i7oPo-|^7ng 
Alternatively, eliminating 6T, we find 

^^^^ = Vf/o ■ 5v„ + '^gn^oV ■ 5v„ - gV ■ (rico^vj 
at 3 

|Po + 2^noneo + f 7o^ngo \ 

WZ^2 lMii2- (85) 

We next make use of the expression for 5ri2 given in (^) to rewrite (|85|) as 

= Vf/o ■ 5v„ + ^gn^oV ■ 5^,. - gV ■ (n,o5vJ - ^ , (86) 
where the relaxation time for the chemical potential difference is defined as 

V ~ ^12 V §70^0 - Ignl ) ~ (TTu 

Although the relaxation time is proportional to the equilibrium collision time T12 defined 
earher, it has quite a different magnitude and temperature dependence, especially in the 
neighbourhood of Tbec, as a result of the factors rico and a which appear in its definition 
(see Fig. 1 in Ref. [0). 

The results in Eqs. (^6]) and (^) were quoted in Ref. ||I2[ for the special case of a uniform 
gas. Here we can see very clearly that plays a key role in describing the equilibration of the 
condensate and non-condensate chemical potentials brought about by collisions between the 
two components. For a uniform gas with spatially homogeneous fluctuations, the solution 
of (|6D is 

W(i) = 5/idifr(0)e-*/^- , (88) 

which shows that any initial difference in chemical potentials decays on a time scale set by 
r^. Similar considerations will of course also apply to an inhomogeneous gas. Thus, for 
fluctuations which occur on a time scale longer than the relaxation time r^, we expect that 
(5/idiff can effectively be set to zero (i.e., p, = fj,c)- It then follows from (^) that we recover 
the equilibrium relationships among all the thermodynamic variables we have defined. As a 
result, our equations in the ^ limit reproduce precisely the hydrodynamic description 
as provided by the Landau two-fluid equations. 

To demonstrate this reduction explicitly it is again convenient to consider the sound 
modes in a uniform Bose gas (see also Ref. |jl2|]). We proceed as in the analysis of the ZGN 
equations in the previous Section, but take into account the fact that according to (PB[), 
fluctuations in Sfidis are coupled to those of the velocities Svc and 5v„. Introducing velocity 
potentials according to = V0c and 5v„ = V0„, the equations in ( p]) and ( ^2|) can be 
combined to yield 

'^~gf^ ^ fi-n^oVVc + 2ghoV'^(f)n + —Sfi^is (89a) 

= hF^ + '^ano VVn + 2(7n,oVVc - — W • (89b) 
oV \3no J 3r^ 
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These equations in conjunction with (|86|) constitute a complete set of equations for the 
normal modes of the system. We look for solutions having a plane wave form 0c,n(r,t) = 
0c,ne'^''''"~'^*'' . In this case, ( p6| ) reduces to 



Substituting this result into Eq. 



1 — iujTf^ \"^" 3 
, we obtain the two coupled equations 



(90) 



gnco 1 



tUJT, 



I + 2gho 1 + 



3fl 



no 



(91) 



and 



'5Po 
3no 



2(7 



n 



cO 



9(1 



?,2 



+ 2^n 



cO 



3(1 - ia;r^) 



(92) 



where the thermodynamic factor a is defined in (0). In the limit cur^ oo, we recover the 
ZGN equations given in (|60D. We also see from (pOD that ^ ignco{4>c — 20„/3)(A;^/c<j), 
which is proportional to k. In the opposite limit ujt^ — » 0, we obtain the equations 



muj~0c = gricoil - cr)k cpc + 2gnQ [ 1 + <; 



3fin 



2gnQ 



I A;^0„ + 25(nco 1 



3rin / 



(93a) 
(93b) 



with (5/idifr — ^ gncoT^{(pc — 2(j)n/S)k^. We see from this that S^^is vanishes in the ut^ 
limit while 6T12, which is proportional to remains finite even when 0. Thus, 

the 5ri2 terms in equations (|5lD and (|52D are still present in the Landau limit, reflecting 
the fact that the establishment of complete local equilibrium requires the continued local 
readjustment of the number of atoms in the condensate and non-condensate. In other words, 
the strong coUisional coupling between the two components, which ensures 5/idifr = 0, implies 
that the condensate and non-condensate densities are not separately conserved as they are 
in the ZGN (cjt^j 00) limit. 

Although not immediately apparent, the equations in ( p3D yield first and second sound 
velocities which are in precise agreement with those determined by the Landau two-fluid 
equations. The latter are usually expressed in terms of pressure and temperature fluctuations 
and yield velocities, u, as solutions of an equation of the form — Au^ + B = 0, where 
A and B are expressed in terms of equilibrium thermodynamic quantities |2^. If these 
thermodynamic quantities are evaluated for our model of the equilibrium properties, one 
can show that the Landau A and B coefficients correspond precisely to those following from 
(^). These differ from the ZGN A and B coefficients given in (|6lD , and as a result, the 
first and second sound speeds found in these two limits also differ. However, as shown in 
Fig. 1, these differences are very small in the case of a weakly interacting Bose gas. That is, 
the transfer of atoms required to equilibrate the condensate and non-condensate is playing 
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a relatively minor role when it comes to determining the magnitude of the first and second 
sound speeds. 

The damping of these sound modes is another matter. The modes are undamped in 
either the ZGN {ujTf^ oo) or Landau {ujt^ 0) limits, but at intermediate values of wt^, 
the damping arising from the equilibration of the condensate and non-condensate can be 
quite significant. To obtain the complex mode frequencies, uj = Q — iT, we can solve either 
the hydrodynamic equations in ( ^T] ) and ( ^2] ) [from (|5BD, we recall that 5ri2 is proportional 
to Sfidis], or Eqs. (|^) together with the relaxation equation ( ^61) for the chemical poten- 
tial difference. These equations yield five mode frequencies: four modes corresponding to 
damped first and second sound {Q = ±Ui^2k) and a purely imaginary relaxational mode. 
The latter is associated with the zero-frequency mode mentioned at the end of Section |V 
and will be discussed further at the end of this Section. 

In Fig. 2 we show the damping of the second sound mode for a uniform gas as a function of 
fir^, for T/Tbec = 0.9 and gn/ksTBEC = 0.2. The relative damping (T/Q) of the first sound 
mode is barely visible on the scale of Fig. 2 and therefore is not displayed. It can be seen 
that the relative damping peaks at Qt^^ ~ 1, which is the transition between the ZGN and 
Landau regimes. The damping of the second sound mode is especially large since it involves 
an out-of-phase oscillation of the condensate and non-condensate density fluctuations. Such 
an out-of-phase oscillation would be expected to have the largest imbalance of the local 
chemical potentials, and thus the largest rate of transfer of atoms between the condensate 
and non-condensate. This damping mechanism will arise whenever the equilibration of the 
condensate and non-condensate is incomplete. 

The results found here for the uniform gas may be relevant to the MIT studies of col- 



lective excitations in highly anisotropic traps [36|, which are believed (on the basis of C22 



collisions) to be in the transition regime between collisionless dynamics and hydrodynam- 
ics. The m = quadrupolar oscillations of both the condensate and thermal cloud were 
observed, as well as an out-of-phase dipole mode. Both the condensate quadrupolar mode 
and the dipole mode involve the condensate and non-condensate oscillating out-of-phase (see 
Section |VII|) . According to the above discussion for a uniform gas, this will entail a large 
chemical potential difference between the two components and we would therefore expect to 
see a large damping rate associated with the equilibration of chemical potentials when the 
number of non-condensate atoms is appreciable. Although such large damping rates were 
indeed observed they cannot be attributed entirely to this source since there will also 



be contributions to the damping from viscosity and thermal conduction which our present 
analysis does not include (see the discussion at the end of this Section and in Section |VII1| ). 
It is nevertheless clear that the relaxation of the local chemical potentials of the two com- 
ponents must be taken into account if one wants to provide a quantitative description of the 
damping processes in condensed Bose gases at finite temperatures. 

In addition to the first and second sound modes discussed above, the ZGN' equations 
also yield a nonpropagating relaxational mode which, in the ^ limit, has a frequency 
uj ~ —i/Tfj., consistent with (|55|). In the opposite limit, 00, collisions between con- 

densate and non-condensate atoms no longer play a role, and the mode becomes the zero- 
frequency mode of the ZGN description. As discussed at the end of Section this mode 
is characterized by the fluctuating variables 6vc = = 0, Suc = —26h and 6P = 2gho6h. 
This mode persists even into the normal phase, in which case the variables Uc and v,, do 
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not exist, and the mode is then characterized by amphtudes satisfying 5P = —2gnQSn and 
= 0. 

In either case, one can show from (^) that the zero-frequency mode is characterized 
by 6P = 0. In other words, even though there is a local change in the total density, 
this mode has a zero pressure fluctuation, and as a result, the fluid remains in mechanical 
equilibrium. This is possible since the change in kinetic pressure, 6P, is exactly balanced 
by a corresponding change in the mean-field potential. It is important to point out that the 
same situation also arises in a conventional fluid having some arbitrary equation of state 
P{n,T). Variations in the pressure are then given by SP = {dP/dn)TSn + {dP/dT)n^T, 
and we see that 6P can be zero if 6n/6T = —{dP/dT)n/{dP/dn)T- Our present model is 
a particular realization of this possibility, and we find from ([7^) that the fluctuations in 
the local temperature and non-condensate chemical potential are related by Sq6T = —fiQdjl. 
The zero-frequency mode in the ZGN description therefore has associated with it a static 
temperature fluctuation. This is consistent from a thermodynamic point of view, but of 
course is not expected in a more complete description of the fluid dynamics. 

Extending the ZGN hydrodynamic equations to include heat conduction will allow the 
temperature fluctuation to relax in time, and as a result, the zero-frequency ZGN mode 
will become the usual thermal diffusion mode with a purely imaginary frequency. However, 
even in the absence of heat conduction, we have found that the zero-frequency mode is 
converted into a relaxational mode within the ZGN' description. It is therefore clear that 
the further inclusion of heat conduction within the ZGN' description will simply combine 
with the collisional relaxation effects due to C12 to produce a hybrid relaxational mode. 
These two effects will also contribute to the damping of the sound modes. In conventional 
discussions based on the Landau two-fluid equations, the thermal diffusion mode found 
above the transition is replaced by damped second sound modes below the transition |7|,|57|. 
Since there is one additional dynamical equation in both the ZGN and ZGN' descriptions 
as compared to the Landau theory, there must necessarily be an additional mode. This is 
the relaxational mode we have been discussing. However, this additional mode is strongly 
damped in the ujt^ — > limit and we then recover the usual results based on the Landau two- 
fluid equations. A more complete discussion of the combined effects of collisional relaxation 
and thermal conduction will be given elsewhere. 



VI. CENTER OF MASS (KOHN) MODE 

In this section, we discuss the center of mass motion of the coupled condensate and 
non-condensate systems. For confinement in an arbitrary harmonic trap potential of the 
form 

[/e.i(r) = ]^m{ujy + uy + ujIz") , (94) 

the center of mass degrees of freedom separate from all other internal degrees of freedom, with 
the general consequence that there exists a special normal mode in which the equilibrium 
density profile oscillates rigidly without distortion. In the electron gas literature, this is 
referred to as the generalized Kohn theorem but the result is independent of statistics 
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and applies to any harmonically confined system at any temperature. The experimental 
observation of this mode provides a convenient way of determining the trap frequencies. 

In the context of a Bose-condensed trapped gas, the theorem implies that the mode 
densities have the form 

nc{r,t) = n^o{r - r]{t)) , (95a) 



n(r,t) = no(r - r}(t)) 



(95b) 



where the center of mass displacement 77 (t) has a time-dependence given by 



(96) 



An approximate theory of the dynamics should, at the very least, be consistent with a mode 
of this type. This is indeed the case for the zero-temperature dynamics of a Bose condensate 
based on the nonlinear GP equation [^, including the TF limit [Q. However, the simplest 
extension to finite temperatures as considered by HZG fails to meet this criterion. In 
fact, it was the small deviation of the center of mass frequency from the trap frequency 
as a function of temperature that emphasized the limitation of the static HFP theory. As 
discussed in Section this approximation does not include the collective dynamics of the 
non-condensate. 



We now prove that our set of coupled equations (^3)^(111) consistent with the gen- 
eralized Kohn theorem. We showed previously [Q that this theorem was satisfied for the 
linearized ZGN equations in (|5TD and (^2|) , which are a special case of the more general 
equations being considered here. However, it is important to demonstrate that this conclu- 
sion remains valid even when collisions between the atoms are included. By doing so, we 
shall be able to clearly see how the effects of interactions, both in the form of mean fields 
and collisions, cancel out in the determination of the center of mass mode frequency. 

We start by considering the condensate density. From ( ^5a| ), we obtain 



dnc{r,t) 
di 



-flit) ■ Vn,o(r - ri{t)) = -V ■ (r7(t)nc(r, t)) . 



(97) 



Comparing this with ( p7a| ), we see that (|95al) is a solution of the continuity equation if the 
velocity field is given by 



(98) 



and provided the source term ri2[/] is zero. This latter result will be demonstrated shortly 
once the form of the distribution function / corresponding to (|95b|) is specified. Using (|98f ] 
in (|27b|), we have 



m 



d^rjait) <9/ic(r,t) 



From the dependence of fic{r,t) on the densities given in (^Sf), we see that (|95|) implies 



(99) 
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/ic(r, t) = jicoiv - r]{t)) + Uextijc) - Uextijc - -nit)) 

= l^cQ + ^rnujlxo.'na,- - ^mujlril, (100) 



since as indicated by (0), /ico is a constant, independent of position. Inserting ( |100| ) into 



(p9|) leads immediately to (0). The condensate equations are thus seen to admit a solution 
consistent with the assumed form of the densities in (P^l). 

A similar analysis of the non-condensate continuity equation (|33a|) shows that ( P5b|) 
implies a non-condensate velocity 

v.(M) = ^. (101) 

again assuming that ri2[/] = 0. This velocity field and the non-condensate density in ( |95b|) 
are clearly generated by the distribution function 

^(p,r,t) = ^^^[_i.(p_^^(,))2+[;„(r_^(i))_^„] _ ^ , (102) 



where [3q and yUo are equilibrium parameters. By noting that ( pSf ) implies 

?7o(r - r){t)) = f/(r, t) + U.^tir - r,{t)) - U.^tir) , (103) 

we see that ( |102| ) is of the local equilibrium form in (|37|), with the non-condensate chemical 
potential given by 

/i(r, t) = ^0 + f/ext(r) - Uext{r - ri{t)) . (104) 

Since jlo = yUco in static equilibrium [see the discussion after (|50D ], this local chemical po- 
tential of the non-condensate is equal to the condensate chemical potential in ( |100| ). 

We now show that ( |102| ) is indeed a solution of our kinetic equation in (|29|). The C22 



collision integral on the RHS of (|29|) is seen to vanish as an immediate consequence of (|102[) , 
namely, that is a local equilibrium Bose distribution. In addition, from (^O]) we see that 
the C12 collision integral vanishes since v„ = Vc and /i(r, t) = fic{r, t). This, of course, implies 
that = 0, as we assumed in the argument leading to (|98|) and (|101|) . A straightforward 
calculation of the LHS of the kinetic equation in ( ^9]) using in (|102|) shows that it vanishes 
if 

= -^\U{r, t) - Uoir - rj)] . (105) 

Using the expression in ( |103| ), we see that ( |105D reduces to the equation of motion in (PB|). 

In summary, we have proven in full generality that the set of dynamical equations in 
(p7|)- (|29|) are consistent with the generalized Kohn theorem for a harmonic trap described 
by (0). In other words, these equations admit a solution consisting of a rigid in-phase 
oscillation of the equilibrium condensate and non-condensate density profiles, the center 
of mass displacement 77 (t) having the trap frequency for a displacement in the a-th 
direction. This result is an important check of the internal consistency of our treatment of 
the dynamics. We also find that this particular solution is, in fact, independent of collisions 



and is therefore valid in the coUisionless |27] as well as in the collision-dominated regime. 



Finally, we note that no assumptions have been made regarding the magnitude of the center 
of mass displacement. These results are therefore valid for arbitrary amplitude of the center 
of mass motion. 
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VII. VARIATIONAL SOLUTION OF THE HYDRODYNAMIC EQUATIONS 



Our purpose in this Section is to develop a practical means of solving hydrodynamic 
equations of the kind we have derived in this paper. For simplicity, we shall restrict ourselves 
to the ZGN equations given in (^Tj) and (|5^), although a similar approach could be developed 
for the more general ZGN' equations. From a practical point of view, the two sets of 
equations give fairly similar results for the frequencies of the first and second sound modes 
in a uniform system, even though they describe different regimes. We would expect a similar 
correspondence in the case of a trapped gas. However, it should be noted that the ZGN 
equations neglect the damping associated with the relaxation of the chemical potentials, 
which is included in the more general ZGN' equations of Section 0. 

A direct solution of the ZGN equations is exceedingly difficult, partly for the number 
of equations involved and partly for the fact that one is dealing with a strongly inhomoge- 
neous system. For this reason it is advantageous to reformulate the problem in terms of a 
variational principle. In doing so, accurate estimates of mode frequencies can be obtained 
by using simple trial functions for the velocity fields whose choice is guided by physical 
considerations. Our immediate objective is therefore to transcribe the ZGN hydrodynamic 
equations into a variational form. 

Although not essential, it is convenient to replace the velocities in terms of displacement 
fields. Such a description is particularly useful for small amplitude oscillations about equi- 
librium, since in this case each fiuid element makes small excursions from its equilibrium 
position and the motion of the fiuid is analogous to that of an elastic medium. The dis- 
placement field for the non-condensate, u(r, t), and the condensate, w(r, t), are defined by 
the relations 

v„(M) = ^. v.(M) = ^?^. (106) 

Assuming a harmonic time-dependence with frequency cj, these definitions are equivalent to 

v„(r) = -iuju{r) , Vc(r) = -iujw{r) , (107) 
and in terms of these variables, the continuity equations take the form 

5rac(r) = -V ■ (ricow) 

Sh{r) = -V ■ (nou) . (108) 

Introducing these definitions into (^) and eliminating the variables Sric, Sh and 6P, we 
obtain 

2 d f dpQ 5 ~ duj\ dUod{hoUj) 



- mriQUJ Ui = —— Uj— h -Pq^ — + 



dxi \ dxj 3 dxj ) dxi dxj 

+ '^d'lT'O g^^^ (noUj + UcoWj) . (109) 

This equation can be simplified considerably by collecting together terms involving similar 
derivatives of Ui, and by making use of the equilibrium condition VPq + ^oVf/o = 0. We 
find 
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d ( du \ d I du ■ \ 

OXi \ OXj J OXj \ OXi J OXiOXj 



X = '^Po + gnl, ^ = Po + ghl. (Ill) 



where 



The first two terms on the right hand side of (|110| ) are reminiscent of the divergence of the 
stress tensor of an anisotropic elastic medium. To make this correspondence more apparent, 
we introduce the symmetric and antisymmetric strain tensors 



1 / dui duj 



Uij = - I h 



2 \dxj dxi J 
1 ( dui duj\ 

• (112) 



2 \dxj dx 

The first two terms on the right hand side of (|1 10|) thus become 



dxi dxj dxj 

The part of this expression depending on the symmetric strain tensor will be recognized as 
representing an elastic medium, with A and fi playing the role of Lame constants . 



We are now in a position to define energy functionals for the non-condensate. First, we 
introduce the kinetic energy functional 

K^[u] = ^Jd^r mfiQU^ . (113) 

Its variation, 5Kn/6ui = mfioUi, gives the coefficient of u"^ on the left hand side of ( |110[ ). 
Similarly, one can easily check that the variation of the potential energy functional 

Un[u] = \j(^^r + l^^l - l^^l + UiUjho ^^^^ (Uext + 2gn^o) | (114) 

yields the desired Wj-dependent terms on the right hand side of ( |11(]| ). The two terms 
proportional to n can be written in several different ways: 

f^i^l - u^^j) = fiul - ^^(V X u) ■ (V X u) 



2' 

dui duj 



dxj dxi 
2 9u, dui 



The first equality in particular shows that the antisymmetric strain tensor is associated with 
the solenoidal part of Uj. 
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The last term in (|110|) accounts for the force exerted on the non-condensate by fluctua- 



tions of the condensate. It can be represented by the interaction energy functional 

Ucn[uM=29jdr^^ (116) 

We have here indicated explicitly that Ucn depends on both of the displacement fields. 

We next consider the condensate. With the elimination of 6nc and 6n from Eqs. (p2D, 
we arrive at the equation 

d d d"^ d"^ ^ 

Here we have defined the operator 



with Cq{y) as given in (0). We emphasize that (|117|) is equivalent to the linearized form 



of the time-dependent GP equation in (||) (of course with the last two terms excluded). In 
particular, we do not make the TF approximation which neglects the quantum mechanical 
kinetic energy. 

The energy functionals for the condensate analogous to those of the non-condensate are 

K,[w] = ^jd^rmnaw^ (119) 



and 



TT \ 1 1 f j3 d{ncoWi) y d(ncoWj) 



The last term in (|117|) follows from the variation of ( |116| ) with respect to Wi. 



We now define the functional 



f/[u, wl 



where 



f/[u, W] = Un[vi] + f/c[w] + f/,„[u, W] 

ir[u,w] = ir„[u]+/i,[w]. (122) 

The requirement that the functional J be stationary with respect to variations of Ui and Wi 
leads to the equations 

5U _^2^K _ ^ 

f-^'f-O, (123) 
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where we have identified cu^ with the stationary value of the functional J. These equations 
are identical to (|110|) and ( 117 ), which confirms that variation of the functional J leads to 



the desired hydrodynamic equations. More importantly, the functional J offers a convenient 
variational method of estimating the normal mode frequencies of a trapped Bose gas. 

One useful approach is the Rayleigh-Ritz method, in which the velocity components are 
expanded in some complete set of functions. This method leads to expressions for the kinetic 
and potential energy functionals which are quadratic forms in the expansion coefficients. 
Variation of J with respect to these coefficients results in a set of linear equations which 
schematically have the form 

\U^P - u'^K^p] = , (124) 



where Uap and Kap are matrix elements of the energy functionals. A simplified version of 
this method will be used to determined several of the normal modes of interest. 



We now establish some general properties of the solutions to Eqs. (|110|) and (|117|) . For 



this purpose, it is convenient to define the six-component displacement vector 

"=(:)■ (126) 

Eqs. (11101) and ( |117| ) can then be combined into one matrix equation 

Lr] = u^Dr], (126) 
where the matrix L has the block structure 

L=(f' f'). (127) 

The elements of this matrix are three-by-three matrix operators defined as 

d I du ■ \ 




dxj \ dxi J dxidxj 



OXiOXj 



(L22w)j = -rico-^iK + g) — {ncQWj) . (128) 



d . d 

Similarly, the matrix D is block-diagonal (-D12 = D21 = 0) with elements 



Dii = mnol 

D22 = mn,oi . (129) 

Here, 1 is a three-by-three unit matrix. From the definition of L, it can easily be checked 
that 
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(130) 



which shows that X is a Hermitian operator. This property can be used in ( |126| ) to show 
that the eigenvalues a;^ of L are real and that the corresponding eigenvectors rja satisfy the 
orthonormality relation 

ld'rril-{Drip) = 5^p. (131) 



A. In-phase and out-of-phase dipole oscillations and the center of mass mode 

As our first application of the above formalism, we consider the center of mass motion of 
a gas trapped in a parabolic potential given by (p^. According to our earlier discussion in 
Section VI, we expect a solution whereby the total density oscillates harmonically without 
distortion. One can easily show that a uniform displacement Ui = Wi = A of both compo- 
nents along one of the Cartesian directions is a solution of ( |110| ) and (|117|) , with frequency 



Lu = LUi, i = X, y, z. According to ( [108| ), this implies a density fluctuation 5nc(r) = —A ^^^ 

for the condensate and bh{j) = —A^^ for the non-condensate. These forms follow from 

(p5|) on expanding in the displacement rj. 

This center of mass mode has the condensate and non-condensate oscillating in phase 
with the same amplitude. One would also expect a dipole mode in which the two components 
oscillate out of phase. In this case, there is no reason to expect the mode densities to be 
given by the equilibrium profiles as for the center of mass mode. However, in obtaining a 
variational estimate of the mode frequencies which is accurate to second order in the error of 
the mode density, it is reasonable to assume constant displacements for the two components, 
but in general with different amplitudes. Thus we take Ui = An and Wi = Ac for a particular 
direction i. For these displacements, the equations in ( |123| ) become 



dU odK 



OJ ^ = 



dAc dAc 

Using the constant displacements to evaluate the energy functionals, we obtain 
K = ^McAl + ^MnAl 

U = \McUjIAI + ^MnUJ^^Al + h,{Ac - Anf , (133) 



where Mc = mNc-, Mn = mN and 



,3 / dnco\ I dn(] 
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(Here, the repeated index is not summed.) This energy corresponds to two one-dimensional 
oscillators confined in an external parabolic potential of frequency uji and coupled together 
by a spring of force constant fcj. 

Substituting (|133|) into ( p.32|) , we obtain the pair of coupled equations 

where = ki/Mc and UJ2 = ki/Mn- One solution of ( |135|) is u = uJi with Ac = An, which 
is the center of mass mode discussed more generally in the previous section. However, the 
present analysis gives another solution 

2 2 I 2 I 2 

U = + t^i + ^2 

with the corresponding amplitudes satisfying McAc + MnAn = 0. This implies that this 
mode is essentially an out-of-phase oscillation of the condensate and non-condensate and 
that the center of mass remains stationary. To the extent that the displacement fields for 
this mode are indeed uniform, one finds that Nc/N = —An/ Ac and hence a measurement of 
the relative amplitudes would provide a direct estimate of the relative number of atoms in 
the condensate and non-condensate. 

These results for the in-phase and out-of-phase oscillations of the condensate and non- 
condensate were previously reported in ZGN |Q. The first observation of such an out-of- 
phase dipole oscillation was made by the MIT group for an anisotropic cigar-shaped 
trap. They find that this mode has a frequency slightly below the axial frequency of the 
trap, and is strongly damped at higher temperatures. Our prediction would give a frequency 
slightly above the trap frequency since the slopes of the condensate and non-condensate tend 
to have opposite signs (see Fig. 3) in the region of space contributing to the integral in ([1361) . 
It is possible that the damping of the mode not included in the ZGN equations would have 
to be taken into account in order to resolve this discrepancy. 



B. Monopole and Quadrupole Modes in Isotropic Traps 

We next consider some of the other low-lying modes of the trapped gas. For simplicity, 
we shall work with an isotropic trap with a common oscillator frequency ujq in all directions. 
The particular modes we shall consider here are the monopole and quadrupole modes, which 
respectively have / = and / = 2 angular character. The lowest monopole mode is expected 
to be a breathing mode whose mode density would have the approximate form A~^no(r/A) 
where A is a time-dependent, but spatially independent, scaling parameter. The equilibrium 
situation corresponds to A = 1 and small deviations from this value imply a displacement 
field of the form u oc r. The quadrupole mode can be viewed as an anisotropic scaling solu- 
tion with / = 2 symmetry. This requires a displacement field with components {x,y, —2z). 
Both the I = and I = 2 modes can be represented by displacement fields of the kind 

u = (ax, ay, hz) (137a) 
^sf = (cx, cy, dz) , (137b) 
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where a, b, c and d are constant variational parameters. For the monopole mode, we have 
a = b and c = d, while for the quadrupole mode, b = —2a and d = —2c. The various terms 
required in the evaluation of the energy functional J for these modes are given in Appendix C, 
where the expressions are reduced to the form of radial integrals over equilibrium quantities. 

One advantage of a variational approach is that the calculated mode frequencies can 
be improved systematically by introducing trial displacements which are more general than 
those given, for example, by ( [1371 ). In the case of the monopole modes, displacement fields 
of the form u = r J2^^q a^r'^'^ were examined |^T|, where the number of terms n in the 
expansion was varied from 1 to 5. The expansion coefficients (together with a similar 
set for the condensate) play the role of the variational parameters, and the functional J can 
be minimized with respect to this set. The total number of collective modes generated in 



this way is 2n, with the lowest pair corresponding to the modes found using ( |137| ). One can 
examine the convergence of the mode frequencies with n and one finds that the results for 
the lowest monopole modes are already quite good with n = 1. We therefore have confidence 
in the accuracy of the present results obtained using the simplest variational form given in 



In Fig. 4, we present our results for the lowest / = and / = 2 modes in an isotropic 
trap with radial frequency = 200 Hz. As for the case of the dipole modes discussed 
earlier, the restricted form of the variational displacement fields gives rise to two modes for 
each symmetry. The calculations were performed for 5000 ^^Rb atoms having an s-wave 
scattering length ofa=110aB. In the lower part of Fig. 4 we show the condensate fraction 
as a function of temperature. The transition temperature is found to be Tbec ~ 149 nK. In 
Fig. 3, we show the condensate and non-condensate radial densities at T = 100 nK, where 
the number of condensate and non-condensate atoms is approximately equal. The depletion 
in the non-condensate density at the centre of the trap due to the repulsive interaction with 
the condensate in the overlap region is clearly evident. These results are very similar to those 
obtained within the static HFP approximation [0] , and confirmed by other calculations 



In the upper part of Fig. 4 we show the mode frequencies as a function of tempera- 
ture. The results below Tbec can be viewed in a first approximation as a superposition of 
the condensate modes given by the HFP scheme [Q] and the non-condensate modes found 
above Tbec EH' extended to lower temperatures. At T = 0, the condensate modes start at 
oj ^ 2.250^0 for / = and cu ^ I.Scuq for / = 2, which are close to the TF limits {^/Eujq and 
\/2ujq, respectively p2|). Apart from some hybridization effects to be discussed, these two 
modes basically follow the HFP behaviour with increasing temperature, merging together 
as T approaches Tbec- The two non- condensate modes are analogous to the / = mode 
at u; = 2uJo and the / = 2 mode at to = V2uJo, as found by Griffin et al. above Tbec- 
However, unlike this earlier work, the present results include the effect of interactions, and 
even above Tbec-, the mode frequencies are shifted slightly from their non-interacting val- 
ues, although it turns out that the / = 2 mode frequency (2~) remains at \/2ujq within the 
present variational approximation. For temperatures below Tbec, the non-condensate den- 
sity fluctuations couple with the condensate fluctuations and these modes show an increasing 
frequency shift with decreasing T. In addition, hybridization of the modes is observed at 
points were the uncoupled modes cross. For example, the u = 2ujq mode is seen in Fig. 4 to 
hybridize with the condensate mode just below Tbec, and once again at lower temperatures. 

This hybridization is most clearly revealed by looking at the amplitudes of the condensate 
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and non-condensate fluctuations associated with each of the modes. In Fig. 5 we show the 
amphtudes for the two / = modes, where the upper / = mode is labelled 0"*" and 
the lower, 0~. Starting at low temperatures, the O"*" mode is essentially a non-condensate 
oscillation (|An+| ~ 1) with a very small condensate component. After passing through the 
hybridization point near T ^ 35 nK, the upper / = mode changes over to a condensate 
oscillation, with a small out-of-phase non-condensate amplitude. At the second hybridization 
point near 136 nK, the amplitude of the non-condensate grows dramatically as the mode 
switches over to the 2ujo non-condensate mode. In the temperature range between this 
hybridization point and Tbec — 149 nK, the condensate amplitude is still appreciable; 
however, the condensate has a minimal effect on the non-condensate oscillation since the 
condensate fraction is quite small here. One can simply view the condensate oscillation as 
being driven by the much more massive non-condensate component. The amplitudes for 
the lower frequency / = mode (0~) are shown in the lower part of Fig. 5. Above the 
T ~ 136 nK hybridization point, this mode is seen to be essentially a condensate oscillation 
~ 1), with a very small non-condensate amplitude. A mode of this type is to be 
expected since the non-condensate is relatively massive and the condensate simply oscillates 
in the presence of a static non-condensate density. In this respect, this mode is clearly 
equivalent to the mode found in the HZG approximation 0, which ignored the collective 
dynamics of the non-condensate completely. 

The behaviour of the I = 2 amplitudes shown in Fig. 6 follows a similar pattern. The 
high frequency / = 2 mode (2+) starts off as a non-condensate oscillation at low T and then 
switches over to a condensate oscillation all the way up to Tbec, since in this case there 
is no second hybridization point. The lower / = 2 mode (2~) is mainly a non-condensate 
oscillation for all temperatures, although a small condensate amplitude is evident as well. 
Interestingly, the condensate amplitude goes to zero at T 125 nK, which implies that at 
this temperature, the condensate is effectively stationary in the presence of the oscillating 
non-condensate. A similar zero-crossing can be seen in Fig. 5 in the case of the lower 
/ = mode. This behaviour can perhaps be understood most easily in the case of the 
breathing mode, which corresponds to a dilation of the density. By referring to the non- 
condensate density in Fig. 3, one can see that the corresponding density fluctuation, 6h oc 
3fio + r{dnQ/dr), will be both positive and negative in different regions of space. Thus the 
zero-amplitude condensate oscillation arises when the radial force on the condensate due to 
the breathing non-condensate averages to zero. Although the situation for the / = 2 mode 
is not as easy to interpret, presumably the force on the condensate is being averaged to zero 
in a similar way. We conclude that a simultaneous measurement of the condensate and non- 
condensate amplitudes as a function of temperature should reveal a particular temperature 
at which the condensate amplitude passes through a minimum. This would be an especially 
interesting experimental signature of the modes we have been discussing. 

We have checked that the results shown here for 5000 atoms are in fact qualitatively 
similar to those for other values of A^, from a lower limit of 2000 to an upper limit of 
20,000 atoms. There is no reason why this qualitative behaviour should not persist to 
arbitrarily large A^, where a hydrodynamic description would be expected to be valid. Thus 
the numerical calculations we have performed can be viewed as representative of the modes 
in the hydrodynamic regime. 

We are not claiming that a trapped gas containing 5000 atoms, approximately the number 
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of atoms in the experiments of Jin et al. |^3[, is actually in the hydrodynamic regime. The 
crossover into this regime occurs at ut^ 1, where Tc is the mean time between collisions. If 
one uses the classical gas expression = anv (where a = Sna^ is the s-wave cross-section, 
n is the maximum non-condensate density, and v is a characteristic thermal velocity), one 
obtains for the JILA experiment |^ the estimate ujTc ~ 4, which places this experiment in 
the coUisionless regime. However, this estimate based on the classical collision time is not 
entirely secure. Our recent work [|r^ has shown that collision times defined by the C12 and 
C22 collision integrals can in fact be much smaller than in the region around Tbec- Thus 
there is the possibility that the JILA data near Tbec is actually closer to the hydrodynamic 
regime than would be expected on the basis of the classical gas expression for t^. In any 
case, it would be very useful to have further experimental studies of collective modes at 
finite temperatures in traps containing a much larger number of atoms. Performing such 
experiments in isotropic traps would further facilitate a comparison with our theoretical 
predictions, not only with regard to mode frequencies, but also with respect to the relative 
amplitude of the condensate and non-condensate oscillations. 

Finally, we comment on the results of Shenoy and Ho which were obtained using the 
Landau two- fluid hydrodynamic equations. Above Tbec, they find hydrodynamic modes in 
agreement with our results, namely, an Z = mode at 2ujq, the center of mass mode at loq 
and an Z = 2 mode at \/2uJo. Below Tbec, these modes continue as non-condensate-like 
modes (i.e., the non-condensate has the larger amplitude as indicated in Figs. 5 and 6) 
and their results are again in agreement with ours. However, there are large differences 
in the condensate-like oscillations which, in both approaches, are essentially out-of-phase 
oscillations of the two components. Referring to Fig. 4, we see that our / = and / = 2 
condensate modes are essentially the modes obtained from the T = time- dependent 
GP equation, but with the condensate number given by the finite-temperature value Nc{T). 
These modes then become degenerate as T approaches Tbec, with a common value close to 
2ujq. In the work of Shenoy and Ho, however, the lowest / = condensate mode frequency 
falls steadily from 2a;o at Tbec with decreasing temperature and gives no indication of 
approaching the expected GP frequency of about 2.25c(Jo at low temperatures. Their lowest 
/ = 2 condensate mode starts out at V^luq at Tbec and stays close to this value with 
decreasing temperature. 

As a partial explanation of these differences, we recall that our results are based on the 
ZGN equations which assume only partial local equilibrium and therefore are not in the 
Landau regime studied by Shenoy and Ho. However, as shown in Fig. 1, the results for 
the first and second sound mode frequencies of a uniform gas are in fact very similar in the 
two regimes, which would lead one to expect the differences also to be small in the case 
of a trapped Bose gas. Perhaps the major difference can be accounted for by the different 
equilibrium properties assumed in the two calculations. As discussed in Appendix B, Shenoy 
and Ho based their calculations on the Lee- Yang equation of state |^^, which is only correct 
to first order in the interactions and, moreover, leads to a uniform non-condensate density 
in the region of overlap with the condensate. In contrast, our calculations make use of 
condensate and non-condensate densities determined self-consistently as shown in Fig. 3. 
It would be useful to recalculate the hydrodynamic modes in the Landau limit with these 
equilibrium densities. 
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VIII. CONCLUSIONS 



In this paper, we have derived a set of equations which describe the dynamics of a 
trapped Bose-condensed gas at finite temperatures. These equations consist of a generahzed 
Gross- P it aevskii equation for the condensate order parameter $(r, t) and a semiclassical 
kinetic equation for the excited atom (non-condensate) distribution function /(p, r, t). By 
hmiting ourselves to higher temperatures, we arrive at a simple and intuitive picture in which 
the excited atoms behave as particles moving in a dynamic Hartree-Fock field. Collisions 
between all the atoms are included and in particular, allow for the transfer of atoms between 
the two components. 

Most of our discussion is devoted to the hydrodynamic regime and the derivation of a 
closed set of generalized hydrodynamic equations for the two components. These equations 
are based on the assumption that collisions between excited atoms are sufficiently rapid to 
drive the distribution function /(p, r, t) toward the local equilibrium Bose-Einstein distribu- 
tion /(p,r,t), with the consequence that the C22[/] collision integral for the excited atoms 
vanishes. However, the Ci2[/] collision integral describing collisions between the conden- 
sate and non-condensate atoms remains finite, and as a result, our equations can be used 
to describe the situation in which the condensate is not in diffusive local equilibrium with 
the non-condensate. As noted in Section |, this situation was discussed many years ago 
in the context of superfluid helium near the A point, but using a more phenomenological 
approach . Previous microscopic studies P^JTSj] of the hydrodynamics of Bose super- 

fluids have generally ignored this possibility and simply assumed that the superfluid and 
normal fluid were in local equilibrium with each other. 

We flnd that the equilibration process associated with the transfer of atoms between 
the two components leads to the existence of a new relaxational mode characterized by the 
relaxation time r^. If we consider the dynamics of the system on time scales set by the 
frequency u, we flnd that the linearized form of our theory is equivalent to the Landau 
two-fluid equations in the limit ut^ -C 1. In the opposite limit, ut^ ^ 1, the Cu collisions 
can be neglected and we recover the ZGN limit recently considered by the authors 0. 
The transition between these two limits can be analyzed straightforwardly for the case of a 
uniform Bose gas. As shown in Section 0, we flnd that the flrst and second sound modes are 
coupled to the new relaxational mode and that this coupling gives rise to anomalous sound 
absorption near ut^ ^ 1. Of course, a complete theory of sound absorption must also take 



into account other dissipative mechanisms such as viscosity and thermal conduction |^ 
We briefly discussed the qualitative effect of a flnite thermal conductivity at the end of 
Section 0. 

One outcome of our study has been the realization that the various collision and re- 



laxation times which appear are strongly temperature dependent. In Ref. []12|, these were 
evaluated for a uniform Bose gas, where it was found that ri2 and r22 become very small as 
T Tbeci but that becomes very long. It is therefore possible to have a situation in 
which UT22 ^ 1, as required by our local equilibrium assumption, but with ut^ ^ 1. These 
conditions will always occur close enough to Tbec and deflne the regime in which the ZGN 
equations can be used. The Landau regime will only be valid at intermediate temperatures 
where ~ r22 and cur^ <C 1. In a separate paper, we shall give a detailed discussion of how 
these relaxation times depend on position within a trapped Bose gas. Since the the edge 
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of the condensate is effectively in tlie critical regime, we anticipate anomalous collision and 
relaxation times in this region. 

In Section |V11| , we developed a powerful variational approach to solve our generalized 
hydrodynamic equations. Although we restricted our discussion to the collective modes 
in isotropic traps, the method can also be applied to anisotropic traps, where a direct 
solution of the hydrodynamic equations (as in Ref. would be extremely difficult. Using 
this variational approach to obtain solutions of the ZGN equations, we found significant 
differences in some of the low-lying hydrodynamic oscillations from the results of Shenoy and 
Ho based on the Landau two-fiuid equations. At present, we believe these differences 
arise largely from the use in Ref. of the TF approximation in conjunction with the 
Yang-Lee |Q free energy approximation for the equilibrium thermodynamic properties (see 
Appendix B). To confirm this, however, will require detailed calculations based on the full 
ZGN' equations. Such calculations will have the added benefit of allowing us to explore the 
full range of behaviour from the ZGN limit to the Landau two-fiuid regime. 

There are other possible generalizations and extensions of the work presented in this 
paper. We briefiy discuss some of these, partly for the purpose of pointing out some of the 
limitations of the present work. Perhaps the most obvious limitation of our simple model 
is the use of the Hartree-Fock approximation, as given in (P^, to describe the dynamics 
of the thermally excited atoms. As noted at the end of Section this corresponds to 
thermal excitations having an energy Ep{r,t) = /2m + gnc{r,t), which is the correct 
high momentum limit of the Popov approximation for the excitation spectrum [^. This 
emphasizes the fact that our kinetic equation (|29|) can only be used to describe low energy 
collective modes at finite temperatures, where ( plf ) gives a good description of the dominant 
elementary excitations determining the thermodynamic properties. 

In the region of very low temperatures, the HF spectrum (53) we have used in this paper 
is clearly inadequate. A first step in improving our treatment would be to use the full 
Popov approximation for the excitation spectrum. In this approximation, the excitations 
are identical to the Bogoliubov spectrum, with the exception that the condensate density is 
now temperature dependent. Indeed, KD [13| have used this Popov spectrum in the case of a 
uniform gas to derive a kinetic equation valid at very low temperatures. The main difference 
this leads to in the C12 collision integral is the appearance of Bose coherence factors which 
result from the transformation from a description in terms of excited atoms to one in terms 
of Popov-Bogoliubov quasiparticles. 

In this connection, we note that generalized kinetic equations can be derived using the 
well-known Kadanoff-Baym method [0], which starts from a specific self-energy approxima- 
tion for non-equilibrium single-particle Green's functions. The full, self-consistent Hartree- 
Fock-Bogoliubov (HFB) approximation for the thermal excitations includes the off-diagonal 
self-energy of the non-condensate (m) neglected in the Popov approximation used here. 
(See Section 6 of Ref. |^ for a recent attempt to classify the various approximations used 
in the treatment of static and dynamic properties of Bose-condensed gases.) As is well- 
known in the Bose gas literature [^,^], the full first order HFB does not give a consistent 
"gapless" approximation of the thermal excitations, since the off-diagonal self-energies are 
second order in the interaction g. Nevertheless, the kinetic equation based on the HFB 
spectrum does allow one to discuss collective modes which are consistent with conservation 
laws. The HFB approximation was recently used to derive coupled kinetic equations 
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for a Bose-condensed gas in the collisionless regime. Work to include collisions within this 
scheme is currently in progress. 

Turning to another basic limitation of our present analysis, we recall that the derivation 
of the ZGN' hydrodynamic equations depended on the assumption that the excited atoms 
were in local equilibrium with each other, i.e., /(p, r, t) ^ /(p, r, t), where by definition 
C22[f] = 0. The next step would be to look at the effect of small deviations from this local 
Bose distribution, / ~ / + 6f. This procedure parallels the Chapman-Enskog treatment of 



transport processes in classical gases. It was implemented recently for trapped Bose gases |30 



to generate corrections to the ZGN hydrodynamic equations due to shear viscosity and 
thermal conductivity. This gives explicit expressions for hydrodynamic damping associated 
with the Uehling-Uhlenbeck collision integral C22 in ( p3a| ). Similar corrections to the ZGN' 
equations will be discussed elsewhere. 

Finally, it should be noted that the situation in a trapped Bose gas is quite different 
from that in a uniform gas. Due to the decreasing density in the tail of a trapped thermal 
cloud, the assumption of local equilibrium enforced by collisions must always break down. 
This leads to the interesting result that the damping of certain collective modes above Tbec 



is governed by the low-density tail of the thermal cloud. It was shown in recent work |]31|j30 



that a good estimate of this damping could be obtained by introducing a spatial cutoff in the 
hydrodynamic equations. Although this type of calculation can also be extended to below 
Tbec, it is clear that a more rigorous solution of the kinetic equation in (^Qf ) is needed in 
order to deal with the kinetics in the low-density tail of a trapped Bose gas. 
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APPENDIX A: DERIVATION OF COLLISION INTEGRALS 



In this Appendix, we give a more detailed derivation of the expressions for the collision 
integrals given by Eqs. ( |23a|) and ( P3b|) , as well as the three-field correlation function given 
by (|25|) . We closely follow the approach of Kirkpatrick and Dorfman |jl3| in their calculation 
of Ci2 and C22, apart from the fact that we work in the lab frame as opposed to the local rest 
frame. In order to obtain a closed kinetic equation for the distribution function f{p,r,t), 
we must express higher order correlation functions in terms of /(p, r, t). For this purpose, 
we treat H'{t) as a perturbation to the zeroth-order Hamiltonian Ho{t), as defined in (|14|). 
We shall effectively calculate collision integrals to second order in g, while explicitly keeping 
interaction effects in excitation energies and chemical potentials only to first order in g. 

The formal solution to ( [T^ ) can be written as (with ^ = 1) 



p(t,to) = Soit,to)pito)Slit,to)-t / dt'Soit,t')[H'it'),p{t\to)]Slit,t'), (Al) 

J to 

where the unperturbed evolution operator S'o(t, to) satisfies the equation of motion 



^^^^ = Ho{t)So{t,to). (A2) 



The solution of this equation is 



Soit,to) = Texp 



I f'dt'Hoit'] 

Jto 



(A3) 



where T is the time-ordering operator. Iterating ( |Al| ) to first order in H' , one has 
p(t,to) ^ So{t,to)p{to)Sl{t,to)-^ f dt'SoM[H\t'),So{t\U)p{h)Sl{^^ (A4) 

Jto 

Using this in (p^, the expectation value of an arbitrary operator 0{t) can be expressed as 
{d)t = (0(t)) = Trp(to)|4^(t,to)0(to)5'o(t,to) 



-I / rft'^i(t',to)[5'o^(t,t')0(to)^o(t,t'),^'(^')]'^o(t',to)}. (A5) 



We first consider the three-field correlation function {ip'^ipip) , which appears in 
Referring to (|A5|) , we can see that there are two ways in which such an anomalous average 
can arise. The first term on the RHS of (|A5|) gives the initial correlation, while the term 
proportional to H' defines those correlations (to lowest order in g) which build up in time 
as a result of the Bose broken symmetry. It is the latter which we identify with collision 
processes. Thus, for the purpose of evaluating higher order correlation functions, we assume 
that the initial density matrix is such that the initial correlations vanish. This is the basic 



assumption we make in the following |jT3|]. Since the evolution according to So{t, to) is number 
conserving, this assumption implies that in evaluating various correlation functions we need 
retain only those terms which preserve the total number of particles. Thus, to lowest order 
in the interaction g, the three-field correlation function is given by 
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(^t(r,t)^(r,t)^(r,t)) 

= -^Trp(to) f rf^'^^(^^^o)^5'^(t,^OV'^r,^o)^(r,^o)^(r,^o)5'o(t,tO,^^(0 + ^3^ 

•Ito L 

xS'o(t',to), (A6) 

where it is understood that only those parts of H[ and H!^ are kept which yield an equal 
number of creation and annihilation operators. 

Let us first consider the contribution coming from H[{t'), which is defined in The 



evaluation of ( [X^ ) and similar correlation functions is facilitated by two key assumptions |]T3 



(i) that the effect of H[(t') in the interval to ^ ^' ^ ^ is essentially a collision process, which 
occurs on a time scale much shorter than all other time scales in the problem, and (ii) 
that the hydrodynamic variables vary slowly in space and time. One therefore expects the 
dominant contribution from the commutator of field operators ?/'(r, t) and terms contained in 
the effective Hamiltonian to come from values of r' and t' close to r and t. In this situation, 
we can expand the quantities nc, n, 9 and U that appear in i^eff about this point. It is 
sufficient to use 



ne(r ~ n,(r,t), n{r\t') - n(r,t), U{r',t') - U{r,t), 
^(r', t') ^ eir, t) + — (f - t) + V0 ■ (r' - r) 

= ^(r, t) - £,(r, t){t' -t)+ mve(r, t) ■ (r' - r) , (A7) 



and 

^o(t,t') -e-'"^"^*^^*"*'^- (A8) 

Introducing the Fourier transform of the non-condensate field operators according to {V 
is the volume of the system) 

^{v, to) = ^7 E ^P^^''^ ^o) = ^7 E (A9) 

W p W p 

and using ( [A7|) and ( |A8|) , one contribution to the commutator in ( |A6| ) is 

[Slit, t')^^(r, toMr, toMr, to)So{t, f), H[{t')] 
2^? 



Pl,P2,P3,P4 



Pl>P2,P3,P4 

^ , (AlO) 



^'^P4,Pc 



'^pl'^P2"P3,P4 + '^pi '^Pa "P2 ,P4 



where we have defined the condensate momentum pc = mvc. Substituting this expression 
into (|A6|) , we obtain 

(V>t(r,t)^(r,t)^(r,t))(i) 

= 2igh{r,t)^Jn^{r,t)e'^^'^'^^ ^ e-i(p.+Pi-P2-P3)-r^^ 



V 

Pl,P2,P3,P4 



P4,Pc 



X r [5p3,p,«ap,)i, + 5p„p,(at^ap3),,] , (All) 

Jta 



I to 

AA 



where 



Following KD |T^, we assume that the initial statistical density matrix p{to) has the form 
appropriate for the HF Hamiltonian HQ{t), in which case 

(«ii"p2)*o - (^Pi.Pa/lPi, r, t) . (A13) 

This identification of the expectation value with the distribution function at r and t is based 
on the picture that the collision processes are short-ranged in both space and time. Finally, 
the remaining time integral in ( |A11| ) is performed by setting to ^ —oo and introducing the 
convergence factor e"*^*-*"* with the result 

/* cit'e^(-^+-'i--'2-.-3+i5)(t-t') ^ ^^(^^ + _ ^2 - £3) + iP ( -) . (AM) 

J -00 \£c-\- £1 — £2 ~ £3 / 

It is useful to note that the approximations made in the above calculations can be 
reproduced simply by replacing H'i(t') by 

H[{t') ^ -2(7r^/2nny2^(5p,pJe-'^e'^^(*'-*)e'P-'-ap + h.c.] , (A15) 

p 

where the variables n, ric, Vc, £c and 6 are all evaluated at r and t. In the subsequent calcu- 
lations it is convenient to make similar replacements for the other parts of the interaction 
Hamiltonian. In particular, we have 

m') ^ ^ E V+P.,2pJe-^2^e^2-(*'-*)e^2p-%,ap, + h.c], (Al6a) 

Pl,P2 

^3(0^^ E <5p.+p„p.+p3[e-^V-(*'-*)e^P-^at,^ap,ap3 + h.c.], (Al6b) 

Pl,P2,P3 

H'iit')^^ E ^Pl+P2,P3+P4«ii«i2«P3«P4 -2^^E4"P- (A16C) 

Pl,P2,P3,P4 P 

One thus ends up with expressions for H'iit') which now depend explicitly on (r, t). These 
forms are of course restricted to the evaluation of correlation functions such as ( [A6| ) involving 
the field operators ip{r,t). 

The contribution of H^{t') to the three-field correlation function in ( |A6| ) can now be 
obtained in a similar way using ( |A16b| ). In the course of the calculation we encounter higher 
order correlation functions of the type {a^a^aa), which we express in terms of the normal 
products (a'l'a) by means of Wick's theorem. We thus find 



PliP2,P3 

^ tn III 



'° p'i,p^,p^ 
x[5p2,p^5p3,p;,(a5,iVi)t' + ^P2,p;,(«piVi)t'(«p;,«P3)t' 
+'^P3,p;,(a{,iVi)t'(«i^«P2)t' - ^PuP'Mp'^(^P2)t'{al'^aps)t' 

+ ^P2,P^(«il"P3)*'(4^Vl)t' + ^P3,pJ,(«il«P2)t'(«p^Vi)t'] ' (A17) 
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where (a^^apj)^/ is given by ( |A12| ). The last two terms in ( |A17| ) have a special form. Isolating 
this pair of terms, we have the sum 



P'l.P2'P3 
P'l.P2.P3 



Pc+p'i.Pa+Pa 



+Pi.P2+P3 



'^P2,P^(«k^P3)t'(4;iVi)t' + '^P3,pJ,(«k"P2)*'^4^Vi) 



P2 



P2,P2 



(^Pi^Pa)*' +^P3,P^(«!.i"P2)*'] (^I's^^pOi' 



= nVY^Sp^^p'^ ['^P2,P^(aii^P3)*' +'^P3,P^,(«ii«P2)t 

where in going from the first to the second line we have interchanged the and Pg summa- 
tion variables, and have used ( A13 ) to obtain the last line. When this expression is inserted 
back into ( [A17] ), we see that the last two terms in this equation exactly cancel the contribu- 
tion from H[ in ( |A11| ). We thus obtain the following explicit expression for the three-field 
correlation function 



(^t(r, t)^(r, t)V^(r, t)) = -r27:gny^e'' — ^ 



Pl,P2,P3 



X 

x5 



5{ec + ei -S2-e^) + -P- 



1 



Vr £c + El- 62- e-i 
Pc+Px,P2+P3[/l(l + /2)(1 + h) - (1 + /l)/2/3 



(A18) 



with fi = f{pi,r,t). We emphasize that in deriving ([A18| ), we have assumed that initial 
values of anomalous correlations, such as {cipiCip2)toy vanish. In addition, we have treated 
the system as locally homogeneous, with the consequence that the local HF single-particle 
energies, ep{r,t) = /2m + U{r,t), appear. We finally note that the three-field correlation 
function in ( |A18| ) is explicitly proportional to ^/n^ and the interaction strength g. Thus it 
vanishes above Tbec well as in a non-interacting Bose gas. 

Following the same procedure, one can also calculate the anomalous pair correlation 
function 



m{r,t) = {^{r,t)^{r,t)) 

= \ E («Pi«P.)te^^''^^^^^'^ 

Pl,P2 



The final result analogous to ( |A18| ) is 

m(r, t) = -ing<l>^^ ^ ^Pi+P2,2pjl + /i + /2] 

Pl,P2 



X 



6{si + £2 - 2ec) + -P- ^ 



vr £1 + £2 - 



(A19) 



(A20) 



We note that if ( |A20| ) is evaluated for a static equilibrium situation (v^ = 0, Ec = fico and fi 
equilibrium Bose distributions), the imaginary part vanishes and the real part simplifies to 

dp 1 + 2/0 



mo(r) 



(27r)3 2E,iv) 



(A21) 
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where Ep{r) = £p(r) — /ico is the local HF excitation energy. This result for the anomalous 
pair correlation is consistent with the well-known HFB expression [^] for a uniform Bose 
gas (in the HF limit). 

The real parts of (|A18|) and (|A20|) contribute to the local condensate chemical potential 
in d^). As is well known (see, for example, Section 5 of Ref. [^] and Ref. |^9[), the anomalous 
pair correlation in ([A21|) is ultraviolet divergent. This divergence is removed in the more 
complete Beliaev approximation [^, which includes all self-energy contributions to second 
order in g. Since we do not include renormalization effects consistently to order g'^ in 
the present treatment, we will drop these contributions to the chemical potential. The 
imaginary parts of ( |A18| ) and ( |A2CI| ) contribute to the source term on the right hand side 
of (|5a]) . However, the contribution from m involves the energy-conserving delta function 
6{ei + 62 — 2ec). Taking the momentum conservation condition in ( [A2(]| ) into account, one 
finds that the argument of the delta function in the TF limit is equal to {P1 — P2 Y /4'm + 2gnc- 
Since this expression is positive definite, the imaginary part of rh makes no contribution to 
the source term in (|5aD. 

Using similar techniques and approximations, we next evaluate the collision terms given 
in the right hand side of (|2T|). The C12 collision integral is defined as the contribution from 
the perturbation in ( |A16b| ), 



Cl2[f] 



-^Trp(t,to)[/(p,r,to),^3W] 
q pi,P2,P3 



P2+P3 



e ^^e^P^ "" 



X 



Pl:P + q, 



/2(Op_q/2'3'P2'3'P3)t " <^p2,p-q/2(OpiOp+q/20p3)t 
-'^P3,p-q/2(apiap2ap+q/2)t " h.c. 



(A22) 



Eq. (|A18|) can be used to identify the three-field correlation function (ap^ap2«p3)t- Inserting 
the result in ( [A22| ), we find 



Culf] = 47r5(Vl^ ^ E "^(^c + 61-62- £3) V+pi,P2+P3 

P1,P2,P3 

X [V,P - 5p2,P - 5p3,p][(l + /l)/2/3 - /l(l + /2)(1 + /s) 



(A23) 



The C22 collision term is the contribution from the H'^ perturbation defined in ( [A16c ). 
We note that the second term in ( |A16c|) is proportional to the number operator which 
commutes with any number-conserving operator. We thus have 

C22[/] = -^^Tp{t,to)[f{P,r,to),H'M 



'2V 



E E - 

q Pl,P2,P3,P4 



3«qr 



Pl+P2,P3+P4 



[5p. 



p+q/2 



p-q/2'^P2'^P3'^P4 



P2, p+q/2 



(a 



p-q/2'^Pl'^P3'^P4/< 



-<^P3,p-q/2(Opi 0-^,2 a'P4'^p+q/2)t " <^p4,p-q/2 (o-pi O-pj '2p3 '2p+q/2)t 



(A24) 



According to 



(«ii^P2^P3"p4)t = TTp{to){S^{t,to)al^al^ap.^apJ{t,to) 



to 



'^Pl '^P2 '^P3 '^P4 



,H',{t')]S{t',to)}. (A25) 
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We find that the first term of order g'^ in ( |A25| ), while finite, makes no contribution to ([A24D . 
Thus the C22[/] collision integral is explicitly of second order in g. Using Wick's theorem to 
evaluate the second term in ( [A25D , we obtain the relevant contribution 



('^pi'^P2'^P3'^P4)i 



27cig 



S{6i +£2- £3- £4)5p^- 



P2,P3+P4 



X /l/2(l + /3)(1 + h) - (1 + /l)(l + /2)/3/4 

Inserting ( [A26| ) into (|A24|) , we finally obtain the expression 



(A26) 



C22[f] 



6{ei + 82- £3- £i)5p^- 



P2,P3+P4 



Pl,P2,P3,P4 



P,Pl 



'^p,P2 ^P,P3 '^P,P4] 



+ /3)(1 + U) - (1 + + /2)/3/4 



(A27) 



Replacing the momentum sum (1/V^) I]p by the integral / d'p/ (27r)^, and the Kronecker 
delta function Vd^^o by the Dirac delta function {2'nf5{jp) in (^T§, ( [A25D and ( [M7D , 
we recover the expressions written down in Eqs. (0), ( |23b| ) and ( P3a| ), respectively. As 



discussed in the body of the paper, there is a close relation between the expressions in ( |A18 ) 
and ( [A23| ). This ensures that the separate continuity equations for the condensate and 
the non-condensate [see (p7a|) and (|47a|) ] are consistent with the requirement that the total 
number of atoms is conserved. (We note that KD |]T3| did not obtain the equivalent of these 
equations since they did not explicitly evaluate {ip^^ ipip) .) 



APPENDIX B: EQUILIBRIUM PROPERTIES 



In this Appendix, we briefly discuss the difference between the static equilibrium func- 
tions we use for a weakly-interacting Bose gas and the well-known results obtained by Lee 
and Yang |Q at finite temperatures (T < Tbec)- The Lee- Yang expressions for the ther- 
modynamic functions were derived as first-order corrections in the interaction g to the ideal 
Bose-condensed gas results. We list them here for convenience: 



P = Pcr + ^g{n^ + nl) , 

3-1 
e = -Per + -g{n'^ - nn^r + 2n 

TS = ^Per + ^fi'(^cr - ^^cr 



The associated free energy density is 



1 



cr/ ' 



(Bl) 



f = e-Ts = -Per + ^g{n^ + 2nncr - nl^ 



(B2) 



Here 



P„ 



:g5/2{z = 1) 
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= -^93/2{z = 1), (B3) 

are the ideal-gas values of the pressure and non-condensate density at Tbec- One can see 
that ( |B1| ) and (p2|) explicitly involve only first order corrections in g. To this order, these 
expressions are entirely consistent with our thermodynamic expressions given by (|65D , (|68D 
and (|69|). For example, we have 

- -f^95/2{z = 1) - (^n,— (73/2(1) + ■ • • 

= Per - ghcr{n - Ucr) . (B4) 

Substituting this into (^) with n = fieri we obtain the expression for the pressure in ( |B 1[ ) . 
However, it is clear that our description of the equilibrium properties (based on the HF 
approximation) goes beyond the Lee- Yang description, which is only correct to first order 
in g. 

The implications of this become apparent when these results are considered in the con- 
text of a trapped Bose-condensed gas. Within our semiclassical HF description, the non- 
condensate density is given by 

n = j^93/2{z) , (B5) 



with the fugacity given by 



e^^'^^-^) . (B6) 



Here, /ic is the condensate eigenvalue as obtained from (^) and U (r) = [/ext(i")+25'n(r) is the 
effective potential as seen by the non-condensate atoms. For (|B5| ) to be well-defined, fic must 
lie below the minimum of the effective potential U{r). Because of the localized condensate, 
the non-condensate atoms experience a localized repulsive interaction 2gnc{r) which has 
the effect of depleting the non-condensate density at the center of the trap where the non- 
condensate density overlaps with the condensate. This is evident in Fig. 3 (see also Refs. 
and 1^), and gives rise to corresponding spatial variations in the other thermodynamic 
functions. 

In contrast, in the Lee- Yang description, the non- condensate takes on the critical value 
ficr wherever the condensate density is finite. This is the approximation used by Shenoy and 
Ho |3^ in their calculation of the normal modes of a trapped Bose gas. Working within the 



Thomas-Fermi approximation, their condensate density is given by 

nc(r) = ^(/i - t/ext(r) - 2ghcr) , (B7) 

which goes to zero at a radius Rtf defined by UextiRTp) = — '^gficr- For r < -Rtf, all 
the thermodynamic coefficients required in the Landau equations were evaluated using the 
Lee- Yang free energy density in 
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Above Tbec, Lee and Yang obtain the non-perturbative result for the free energy 
density (generahzed to include a trap) 

/ = --^5'5/2(^o) - gn^{r) + nno - nf/e^i(r), (B8) 

where the density n(r) is obtained from ( P5| ) with the fugacity given by 

The chemical potential fiQ is chosen to give the desired total number of atoms, A^. The other 
thermodynamic functions are given in this temperature range by 

^- ^ g5/2{zo) + gn^, 



(3A^ 
3 1 



e = 2 /3A3^3/2*^^°'' +gn + nUext{r), 
5 1 

sT = -j^gr,/2izo) -nksT In zq. (BIO) 



These are the results used by Shenoy and Ho for the thermal cloud in those regions 
where no < n-cr, that is, for r > R^f when T < Tbec- We note that these expressions take 
interactions into account in a self-consistent fashion, unlike the situation for r < i?TF where 
the non-condensate density takes on the ideal gas value of n^r- 



APPENDIX C: 

We list here the various integrals required for the evaluation of the energy functional J 



in Section |VI1| . The form of the non-condensate displacement is u(r) = (ax, ay, hz) and that 



of the condensate is w(r) = (cx, cy, dz). For the kinetic energy terms, we have 

KJvi] = -{2a^ + &^)m / d^r^noir) 
6 J 

KJw] = i(2c^ + d'^)m I d^rr^ricoir) . 
6 J 

In obtaining these results, we have made use of the spherical symmetry of the equilibrium 
densities. The various potential energy contributions are the terms appearing in (|1 14|) , (|116|) 
and (|T20D : 



lld'rul = ^i2a + bY Jd'rXir) 

^ J dhii{ul - 4) = i(2a2 + 62) j d^rfi{r) 

1 f &^ 1 ( 

- / d^rUiUjfiQ— Uext = 7'(2ct^ + b^)mujQ / d^rr'^fiQ 



J oXiXj 6 



-Jd rUiUjno-^-—2gnco = —{a - b) g J d rrriQ- 

^ 3 



cO 



dr 
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+ ^{8a^ + 4ab + 3b^)g J d^r r^no 



2g [d^r^^^^^ = ^{8ac + 2ad + 2bc + 3bd)g [ d^rr^^^ 
J oxi oxj 15 J or or 



y,.^SJ^,0^ __ 1 ,,,,,, (rf) A f.^- 

-^4 1' 



^3 

.2 „ /a^ \ 2 
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Figure Captions 



Fig. 1: Squares of the first and second sound velocities (normalized by the first sound velocity 
of the ideal gas at T = Tbec) vs. T/Tbec- The solid hnes are for the Landau two- 
fluid theory while the dashed lines are for the ZGN hydrodynamics. The curves were 
generated for gn/kBTBEC — 0-2. 

Fig. 2: Damping of the second sound mode in a uniform gas due to the relaxation of the 
condensate and non- condensate chemical potentials at a temperature T /Tbec — 0.9. 
The damping, relative to the mode frequency, Q2, is plotted as a function of il2r^. 

Fig. 3: Density of atoms as a function of position in an isotropic parabolic trap containing 
5000 ^^Rb atoms. Parameters defining the trap and atomic interactions are given in 
the body of the paper. The solid line is the total density and the dashed line is the 
non-condensate density. The parameter d is the oscillator length, d = Jh/mujQ. 




Fig. 4: Monopole {I — 0) and quadrupole (/ = 2) mode frequencies, in units of the trap 
frequency cjq, vs. temperature for 5000 ^^Rb atoms in an isotropic parabolic trap. The 
'+' and '-' signs merely distinguish the two modes of a given angular character. The 
lower part of the figure shows the condensate fraction as a function of T. 

Fig. 5: Variational condensate [A^) and non-condensate {An) amplitudes for the monopole 
modes as a function of temperature. Panel (a) is for the 0+ mode and panel (b) 
is for the 0~ mode. The vertical line at T ~ 149 nK denotes the BEG transition 
temperature. 

Fig. 6: As in Fig. 5, but for the quadrupole 2+ and 2~ modes. 
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